{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparing collections (Part One)\n",
    "\n",
    "* Set comparison\n",
    "* Ordered collections\n",
    "* Ranked collections\n",
    "* Collection transformation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import random\n",
    "import collections\n",
    "import operator\n",
    "import time\n",
    "\n",
    "import numpy as np\n",
    "import scipy\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn import datasets\n",
    "import twitter\n",
    "%matplotlib inline\n",
    "\n",
    "import count_min\n",
    "\n",
    "# some matplotlib color-mapping \n",
    "cmap = plt.get_cmap('viridis')\n",
    "c_space = np.linspace(0,99,100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Set comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# create two sets of data scientists\n",
    "a = set(['josh','fiona','scotty','skippy'])\n",
    "b = set(['jeff','whitney','fiona'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jkolb/pyenvs/3.5.1/lib/python3.5/site-packages/IPython/core/formatters.py:92: DeprecationWarning: DisplayFormatter._ipython_display_formatter_default is deprecated: use @default decorator instead.\n",
      "  def _ipython_display_formatter_default(self):\n",
      "/Users/jkolb/pyenvs/3.5.1/lib/python3.5/site-packages/IPython/core/formatters.py:669: DeprecationWarning: PlainTextFormatter._singleton_printers_default is deprecated: use @default decorator instead.\n",
      "  def _singleton_printers_default(self):\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'fiona'}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a.intersection(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "c = set(range(10))\n",
    "d = set(range(8,15))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{8, 9}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c.intersection(d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14}"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c.union(d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# List comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Number sequences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# create an ordered sequence of evenly-spaced integers\n",
    "# create a second ordered seqence of intergers, which differs from the first by only noise\n",
    "a = np.array([[i,i+np.random.normal()] for i in range(100)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x114d83fd0>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAECCAYAAAASDQdFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XecVOXB/v/Pmba9sJ2FpcMNgkhRQURABaVYY4wxGp+Y\nR6OGX4qJaeabRFNNnliiicYWe2LUKBYURVBQVFAEpd4U6cvCsrC9z5zfH7MqSQiSYeHM7Fzv14sX\nM2dmz14zO3vtPfc5c47jui4iIpIcfF4HEBGRo0elLyKSRFT6IiJJRKUvIpJEVPoiIklEpS8ikkQC\nh3pHY8wY4CZr7anGmBHA7UA70AJcZq2tNMZcCXwNaAN+Za2dfSRCi4hIbA5ppG+M+R5wL5DSseg2\nYKa19jTgGeAHxphi4BvAScBU4DfGmGDnRxYRkVgd6vTOBuD8/a5fZK1d0XE5ADQDJwJvWmvbrbW1\nwHpgeKclFRGRw3ZIpW+tfYboVM7H13cBGGPGATOBW4FsoGa/L6sHcjotqYiIHLaYN+QaYy4C7gSm\nW2urgFqixf+xLKD68OKJiEhnOuQNufszxlxKdIPtJGvtx8W+BPilMSYEpAGDgZUHW4/ruq7jOLFE\nEBFJZjEX539d+sYYH/AHYAvwjDHGBRZYa280xtwOvNkR6HprbevB1uU4DpWVdTHEjg+FhVnK7yHl\n904iZ4eukT9Wh1z61totwLiOq/n/4T73A/fHnEZERI4ofThLRCSJqPRFRJKISl9EJImo9EVEkohK\nX0Qkiaj0RUSSiEpfRCSJqPRFRJKISl9EJImo9EVEkohKX0Qkiaj0RUSSiEpfRCSJqPRFRJKISl9E\nJImo9EVEkohKX0Qkiaj0RUSSiEpfRCSJqPRFRJKISl9EJImo9EVEkohKX0Qkiaj0RUSSSOBQ72iM\nGQPcZK091RjTH3gQiAArrbUzO+7zU2AG0AZca619t/Mji4hIrA5ppG+M+R5wL5DSsegW4Hpr7UTA\nZ4w51xgzEphgrR0DXAz86UgEFhGR2B3q9M4G4Pz9ro+21r7RcfklYAowHngFwFq7DfAbY/I7K6iI\niBy+Qyp9a+0zQPt+i5z9LtcBOUAWULPf8vqO5SIinmtra8N1Xa9jeO6Q5/T/RWS/y1nAPqAWyP6X\n5dWftaLCwqwYI8QH5feW8nsnUbK3trbyt4euJj9zGY3N6ZT2+xaFE7+QMPk7W6yl/74xZoK1diEw\nDZgPbAR+a4z5PVAGONbavZ+1osrKuhgjeK+wMEv5PaT83kmU7B8sW8BbC+/ie1daUlJ8QAPPz/st\ndXXTaG72Ol3sDucPVqylfx1wrzEmCKwBnrLWusaYN4C3iU7/zIw5lYjIYZr38h0c2+8Bhg2oIyXl\n00mI3qXVVFZWkpVV6GE67xxy6VtrtwDjOi6vByYd4D4/B37eWeFERA6F67rMf/XPtLWtoq29iDPO\n/AG0PsuAvmF27fax/qMWBvaL7ny4fG0/LhlfRnV1Ag/1D0OsI30Rkbjxysu3cfyIeyjMh/Ub2rjv\nrlcoK4lOP508JpXXFzXxyhsZ5OSN5/jx3yYYDALJWfr6RK6IJDzHfZ/CfNi1O8ymj8J8Z+Y+srLr\nsRvaAAiFChk84pdMO+9mikvKPE7rLZW+iCS85pbo3uHvL29l+pQ0AKZOSaWqtpVf3TGC9OJ7GDFq\nipcR44amd0Qk4Y0Z+33+/o/t1NetZk9VhMICPwC9eqUyqvVz9Os31OOE8UOlLyIJJRwO4/dHS726\nei8LFtxCMNhIn/5XMmDAScyd+2v69ZyLC2zZOZXzPneOt4HjjEpfROLatm0beXf5/VTt2Y0/sJlu\nebXY1QH69SmiYud6rr6qDcdxWLlqPps2/ZYLvvB7du/eDcCocUUep48/Kn0RiVtVVZW88d7XmDR1\nJ3Ofb+Dsc1JZuriF88/3E4lUMHiQi+OEABg2tIm5r74GnElRkcr+P9GGXBGJW+++9yLjJ5cDkJYa\nXVZbE6FnmZ+cHB+VlZ8eEaa93SUcTs5DK/w3NNIXkbjT2trKU8/9hB0Vb1PYL0L3nn5q61za2138\nQYea6jA5uX4iDsx+sZWSknQ2bjiO8z93rdfR455KX0Q8N/f1R9hVOw83HKIw43TsxjlMu/hdxqU4\nzHuhla2bAxR1z+a+O3M55pgBPPpIOca00d6eRnHhV+nXbwJjxhTiOM5nf7Mkp9IXEU+9897LhPP/\nyHEnhHlrbh2tOYvIKGomlJIOwOlnZTJ3VhpnTZrF2ad2jxb7WR6HTmCa0xcRT+3YtZQe/cI01IVJ\nS4VBw0KEIw6tLZ8e+z7g60X37qUayXcCjfRFxBPbtm/kxUU3sK18NUXDIuC4pKZFS33c5EzmzamD\n9hyyUo/ltJOu9zht16HSFxFPzHnr14yasYGRbpBFL9ZAW4g920MMPCZCWrqPopJSBhbcxPBhJ3sd\ntUtR6YuIJ3yp0XMsOY7D+Bm5rJzfh29+71FenvcALW17GTZoCmbgSI9Tdj0qfRE5qpa8/yofbP4r\nmzaVM3SCgz/g0NYawdfaB7/fz/QzrvA6Ypem0heRo2ZH+VY+3PM7Bk9upm9TkDl/ryE/cwB5acfx\nxbN/5HW8pKC9d0TkqKivr+fZFx+mz4hGAFLSfJx+cTe6507gks/dQEpKiscJk4NKX0Q6XSQSYePG\ndWzfvo1IJIK1K7hv9mUEBj/P+g+aPrnfnh1QUmA8TJp8NL0jIp2qra2NO/76bdKHrmX76gactiAt\nTQ1M/XIukMGad2uZ/2QtxfkDKE6ZxPhp072OnFRU+iLSafbt28tDT95Mr3M24g+EqFhfw4nn5PDB\nvMZP7jPkhGzWtvr5ymmPkpqa6mHa5KTSF5HDEolEePn1WdiPltOU/yH7QuX0TSmgub6dzNxoxWQW\nhNi4ooH+x2bQVB/BVz1ahe8RzemLSMwikQi3PHwdO/o9zM60V+k3waXPqBxWLNxHamaAPTtbiERc\n+h+XTXW1j3ceLaN9xZe58qKbvI6etDTSF5H/Wjgc5q4nfoHd+R59J7aRmpmGPxg9hEJ+j3Qi7S6v\nP9jIkOJz2TanmZTMMEMzRzB95iUeJ5eYSt8YEwAeAvoA7cCVQBh4EIgAK621MzsnoojEC9d1mf3a\nU7y18iV6nFNDty3NOE70fLX4fOwtbyavNJXUtFROHHAmF5/9DW8Dy7+JdXpnOuC31p4M/AL4NXAL\ncL21diLgM8ac20kZRSRO3PnEr1nb8x80Fm0mJSNA98E5rH1nHy2NYYZMLGT1iy7Vr46k156vqvDj\nVKzTO+uAgDHGAXKANmCMtfaNjttfAqYAzx5+RBGJB+FwmHL/h5TmBMHv0FjTRnpOkBHnlvHOvdVM\nHvEFbvza+eR1y/c6qhxErKVfD/QF1gL5wNnAKfvdXkf0j4GIdBGO40B7dHJgwCkl2FfL8TWmUpZj\nuO5LP8P0H+pxQjkUsZb+tcAca+2PjTE9gNeB0H63ZwHVh7KiwsLEPpGx8ntL+Y+s6upqVq5dgek/\nmOLiQib0Po/Fa/5BZplDfnoZ10y7kWGDh3sdMybx/twfKbGW/l6iUzoQLfcAsMwYM9FauwCYBsw/\nlBVVVtbFGMF7hYVZyu8h5T+ylq16l8dX/YnQoDZaZvmY1v1SzjzpIgZsGklVzXYGTDqOvLz8uH4M\n/0m8P/ef5XD+YMVa+rcBfzHGLASCwA+BpcB9xpggsAZ4KuZUIuIZ13V5e+mbPLHkXrqfFwICZBbB\nvAWzOHXMmfTvO4ixhaMTujSTWUylb61tAC46wE2TDiuNiHjutr/9mt2DLHVZlXSnxyfLI/52D1NJ\nZ9EnckUEgIXvvcZP7/4eW4o/IKMkFV/IT/X2egCa9rSSuiuPdevX4rruZ6xJ4plKX0R4c+nrzKp+\nhL2Dd+BLidZC75NLqdndzJbHmtjzUgvhEzbySMX/49a/3qDiT2A6DINIEmtsbOTeF+5k+bYl9Lyw\ngPRICuuf30B2STqB1AAZjXkM7HEsNWPfI5gWPUBaQ+465i6Yw8ih4z1OL7HQSF8kif3uiV+ybdRm\nmvKaiIQjOD6H/jP6s+rJHRS/M4xvnPxT0rMyCKZ9Oj5MzQmwr3avh6nlcGikL5JEWlpa2LZtK7WN\nNby+bj7r29bRI1BK97E9WfvCRxQOyCVQn8IXjv8fLpzyJQBOD8zgjtfepGiig+u6VL0WYNrXzqal\nxeMHIzFR6Yt0ca7r8tpb81n90SpWtKymOqeB8L56ek3uQ8u8aHMH04L0nT6Q0Mup/PDSn5Cd/ekH\n6rsX9+Cq8Tcyd/EsHBy+fNalZGdna5fNBKXSF+nifvfX/2Nt3+3sadhB78n9aHtjD71P7w1Adp88\nNr6yicz0DIpbS/jWhd/7p8L/WO+efbmi57VHO7ocASp9kS6somInK7LWk1mYhz8U3YQXSA/RvK+Z\ntLw0cvvnkZqTxvTa6Zw1+RyP08rRoNIX6cLC4TCOL3pyE4J+GnbWUTiqlK1zLDmFmaSF0hjcPIQZ\nF5/tbVA5alT6Il3UI7MfZXnNGmrX7iFUmk7J2N7s/sdGykqKOLHvBcwYexZ+v49u3fK8jipHkUpf\npAt6YcFsXsl6j2D/dPJGGHb/fR1nHXsm0748k149enkdTzyk0hdJUI2Njdz4t5vZ7ttHdiSFmRO+\nzLCBx2A3Wl5c+jLB6ZkAOD6H9IklnFxykgpfVPoiiaKxsZFX336N3MxsTjlxPH+YdQ/rR4Vx/Lns\nfH8z1z77C1KbIWVwNvW5dRTU+gllpwGQUu5SNrrM40cg8UClL5IAqqv38a1Hf82ukZm41a2c8JdF\ntKS24/h9NGzZQ2pOCtnHl1H1pqXb0ELS3AJ2vbaOlEiAnlklXDLk8+Tk5Hr9MCQOqPRFEsDDr/6D\nyjHd8PscSA2xpHkvY8qzCTe10lxeTfG4PgCf7KnjOA4lpxkKlznc/KVfeJhc4o2OvSOSACJEPin0\n9rpGalZvoTJSS8nCBvq05FG/rLzjfn4ad9ZGL2+uZ0LpCZ5llvikkb5IAjh/7Jm89dLt1I/Ip+ad\ndZRMHUaF4xApr+eKjAlkpKbzwrL5FKf3pXhjJvm1BRzXZxhjR4zxOrrEGZW+SBxra2vj9qfvp6K9\nluG+UoIr/MzvXYjjREf9vtJM3reWGy66ltPHTPI2rCQElb5IHPvFY3/gHdOMLxgg3NzMhA1ZZDX4\n+fjEhW44QiYhTzNKYtGcvkgc2xjZgy8YHZv5U0NsYi+XDZpKaEklkRWV9FrSwtfPvtzjlJJINNIX\niSPzFi9k8ebVZPtSyE/Pon7XXuDTXS2zIiHOmTCNGSefQVNTI5mZWd6FlYSk0heJA62trTw97wXu\nq1tGW2k2NfMXkDdpGK0Dcwk/v5ycsmJKWlP5xvSrAPD7/Sp8iYlKX8RjL7wxl7s+eIWKqgqypoyg\ncdVmcscZfKlBUssKiJTmMX1LPt+84kqvo0oXoDl9EQ+Fw2Hu/uAV6kZ0J5x24A2yjuMQDGljrXSO\nmEf6xpgfAucAQeBOYCHwIBABVlprZ3ZGQJGurKmpkYbU6O6XaYPLqF64grT+pdS+upLcGSNxAn6K\n36viS/9zlcdJpauIaaRvjJkInGStHQdMAnoBtwDXW2snAj5jzLmdllKkC3Bdl0gk8k/LMjOz6N+U\nghuOEMzLJqtvTz5fV8rsb9/BVQ2DuWxPb+6+/Cc6bo50mlhH+mcCK40xs4As4PvAFdbaNzpufwmY\nAjx7+BFFEt9jLz/LX9e+Q5vP5aSMUm64/Bts31nOvc88zbH5vem5qYGmgMvo7mO58KKzALio+wUe\np5auKNbSLyA6uj8L6Ac8xz+/a6gD/v3syiJJaPPWLdy99V1ah5USaWnlH4uX8eYNV9EUdOGUYyAc\nYeCHe7n36z8jpLl7OcJiLf0qYI21th1YZ4xpBnrud3sWUH0oKyosTOzdzpTfW/Gev6mpiYXvLqCp\nMAs/0Lh0NVkTh7PnfUvGCSZ6OIWAHzswg+VrlzLj9DO8jnzI4v25/yyJnj9WsZb+m8A3gVuNMaVA\nBjDPGDPRWrsAmAbMP5QVVVbWxRjBe4WFWcrvoXjPX15Rzjf++kc25vgJl1eQMW44vtQgjq/jTXHE\nBX/HSctb2mlrdeL68ewv3p/7z9IV8scqptK31s42xpxijFkCOMA1wGbgPmNMEFgDPBVzKpEE57ou\nf3rpKbYO7U7IcWgL+mh6ZQnZjp+w65I2rD+1Cz4gc+wQnOY2Tt4VZNz5OiKmHHkx77Jprf3hARZP\nij2KSNfwxvtL+P1rs9i8ewcpPaLHsw+WFBJsaOf+KZdx6ytPUO42MbKkP+Oa+lBW0oOJ547/5MiZ\nIkeSPpEr0olc1+Xm15+lfGgP3HSX5g1bSR3Qi0hrGyPb0xnUfxB//vpPgMSfYpDEpNIX6ST3Pft3\nXvxoFR811ZBCKSl9e9K6tZzs11dwwehTuOLqizSaF8+p9EUO07x33uTpBXNZUhQAU0Lbou2EIhEc\nn49QtxwuLBzKVed/yeuYIoBKXyQmTU1N3Pvck7y7diXry3Koa60ivfgYANJHDqVh8QoGBbOYdsxo\nvnbexR6nFfmUSl/kEEQiERoa6snMzCIcDnP1H29iZf8immgko6gvgaZGWrfvItSzGF96GiXde3Dz\n1C8zoE8/r6OL/BOVvsi/cF2XurpaMjOz8Pl8vLl0Cb+d+xx7Aw7+LdvIycpia59ign5/dF97INS7\nJ01r1hPZvJNhPftw0fAJKnyJSyp9kf1sL9/B9x+6h62+CHntLj+cfh63zZtNxYAymtauIzByCLXh\nMG5jIwCBgm40rlxHar8y8tOyuPb4qZx3auJ8qlaSj0pfZD83P/M4G/r2AGDzxk1885F7CacGCVCK\n29pKsFv0kFINS7fj+nz4s7MYsKeJS7MGMPqkEfQu6+VlfJHPpNIX2U9dx6GPW3dW4IRC+MaMpuHd\n98gIhwGHSHMLvtQUMkYPJ33Jh9w4ajLjLh6rA6VJwlDpi3Sort5H71AqHzY20r6nivThQwHIOO44\nGpevYoA/ROCDDbQU5JLu+Lnm/EuZdOI4j1OL/HdU+iLAfc88ySNrVtAY9JO9aSv98/PYULEbp6QI\nXyhIXo+e/GT8mZwwYpTXUUUOi0pfkl5VVRUPr1lBS+8yfEBdSQlnRvyc2a0bz69bhePCOUOGq/Cl\nS1DpS1KLRCKsXbuaxpQQ/o5ljs9HY1uYr5xzAV9BZ6+SrkWlL0mrsbGRmbf9jpV+H00fbSCjqBDH\n5yN1125OG3+61/FEjgiVviStP//j76wuLcHv95Oe143mRW9z5rDjmDH+NCacONbreCJHhEpfkkYk\nEuHNxW/T3t7O0EGD2bRzB073YgD8qamEjOGrk6cyZPAQj5OKHDkqfenyFi5ZzOJVK3lz+ftU9OlP\nc0U5AWbR7g/ga6wnpX9/3EgE09jMgP4DvI4rckSp9KVLe+7Vudy6dCn79u0lZaDBn5ZOZNsmQsOH\nEwKay8tJf28pZ4wYzde/8R2CwaDXkUWOKJW+dClPvvgCr621BB2H/506lbmrVtKaV4BbuRtfMPqp\nWcfv++T+qaWl9EtP54eXX+lVZJGjyvfZdxFJDPMWvcmfPljNh2k5vFVTz+W33sriVSsBSO/bn7oP\nluFGIrit7YTr6wHw1VRzcu++XsYWOao00pcu4cNVK3n4hedp7zUQt72d1qpKskaNoaWqkrqVH5Le\ntx+5+QUcs20LQ48/kZDPYV9bK8cOH8nZp072Or7IUaPSl4T399kvcPfSD9nX1E5qYwOR9jaCOXkA\npOQXEsjIZNCubfz8mpn07tXb47Qi3lLpS8J79v3ltOWVkNGtkNpVywj5HCJtrdCzDACf6zLp+BNU\n+CJoTl+6gEjH/47jkDNsFMeWlfHY967juPq9DKjZw+fzs/nK5z7vaUaReHFYI31jTBHwHjAZCAMP\nEv0dXGmtnXnY6UQO4s3FS/jTMy9Rt6sC1/Xj5BUSrK5i+uhjGTLI8Kfrvu91RJG4E/NI3xgTAP4M\nNHYsugW43lo7EfAZY87thHwiB1ResZPrHnqS5cFs6gaNILy7glOaqrhpxmQunDbD63gicetwRvq/\nB+4CfgQ4wChr7Rsdt70ETAGePbx4IlGtra3c/fjjVDc2c/KxQ9m7by/7coo+GbUEBx1L99wQJ44c\n7WlOkXgXU+kbY74C7LbWzjXGXN+xeP93DXVAzmFmE/nEdf93M8v93XAcePrP95OdGqQ9v5hQacc5\naRvq6D34GG9DiiSAWEf6lwMRY8wU4DjgYaBwv9uzgOpDWVFhYVaMEeKD8h95dXV1rKltxlcYoMZ+\nQObwEyEQoG3zOsIrl9KjpIgzhvTnfy9JvI21ifD8/yeJnB0SP3+sYir9jnl7AIwx84Grgf8zxkyw\n1i4EpgHzD2VdlZV1sUSIC4WFWcp/hOzatYsP16zm2CFDyM/LJxhppw1w/H58gejLNrPPIAIVW3ns\nB98nIyMjbh/LfxLPz/9nSeTs0DXyx6oz99O/DrjXGBME1gBPdeK6JYnMX7SIm2fNpS69G5mzF3Lt\nWady6fgxPLBoKZGmBsLNTfhT0wDoHvKRkZHhcWKRxHHYpW+tPW2/q5MOd32S3Kqr9/HAi6/SlNeT\nALB751Z+9vDjDO7Vk6smHM+Qfv14dsFC1u+toiAjxNev+IrXkUUSij6RK3Hj9gcf4vkPNrJ37x5y\nBvegoXwzKflFhLrlswm4f9Ey7hs1iuuvuQZI/LfoIl5Q6Yvn7IYNPD5rFvN2NhIs7IWvvp6W6irC\nTQ2kl/X55H4NGfm8v2IlPUp7eBdWJMGp9MUT23fs4NaHHmdbxS52t0GzEyQlP7oDWHbfY6gv30Re\n815aaqvwZecDkNFQxchh+syfyOFQ6ctRdft9D7Jw1Sa2b99K9tCTqa7ZRvaAYaS3tVCzaSU5g0fh\nOA4FGWn85prvs2TFChas2UjA5/DFqRPp2aOn1w9BJKGp9OWomfPqPF6weyG9FDLqO5ZGD5fmC6aQ\n2XMQdR8u4syxx3PO1OkMHzqU4UOHcoV3kUW6HJW+HHF1dbX8/PZ7WLpyLcGBYwEItzYDkJrfg5oN\nK8judwyBSDsXT57Ad678qpdxRbo0lb4ccb+5834+bMimPbc3bbu3k17Uk7SCHtSsWUz3kiJ69crn\nhMG59OszgpPHjPE6rkiXptKXI66yrgXHySI9rzt1FZuoX/8uJwwdxIWXXMnJY070Op5IUlHpyxG1\n+N13aavfi5uWh+PzkVncBxPK4rb/912vo4kkJZW+dLq9e/cyb8FCXnv7PdY3ZxOhmLZ1b9KrX3+K\ns9P5wVX/63VEkaSl0pdOtXXbNn7wf/eysylIhAAZBd0ACA4cz+i+Aa69ShtpRbykc+RKp3rk6dnU\npJYBLv5gyifLHcchHIn85y8UkaNCpS+dynWj/6flFFO/cyORcDsAGQ3bOWfKqR4mExHQ9I50ktbW\nVu564FEqdu4isqcRX8EAupUNJXP3ck6bMI5pp11Ovz59vI4pkvRU+vKZwuEwW7duITMzi2XLP6Sy\nah9nTJ5EYUEBz86ew5w3lmPXrSW11xh8/lJaI1sZ0LqB0yedxIXnfQufT28oReKFSl8OqrGxket+\n/Ds27/FRvXs9uWUjcHxB7n7sO+RkZ9EUyCe1W09a/Dmk+6Mvp/T8XqRn1nLR587zOL2I/CsNweSg\n7nvoccobisBJIa1bb4IpWVTvWktO73Hsaw6QlleG4zi4kfA/fV1aSC8tkXik30w5qOaWdhyfnwgu\nOA4Aji+I4/MRSs+lqboCgFBGHnu3LKO1ejsF4W1cc1ninaRcJBloekcOqL6+jlv++AAbPtpOc62f\n1Jze1FZuJDUrHzfShhsJk9mtB9UV6wjX76RPzxJOn346UyaNo7CwCL/f7/VDEJEDUOnLAd3wmz+x\ncXcOjtOftvYtZDev4oJpx1NS1I26+lzWb95BdaPLQFPItVdfSs8eOpuVSCJQ6csBle9uwHGin6bN\nzutNYc4+vvvNqzxOJSKHS3P6ckA5GaFPLruuS3ZG0MM0ItJZVPpyQN+8+gsUpFbgb91KWe4evv31\ny7yOJCKdQNM78k/+9vdnmDN/Oa4LE8cN5qv/80Wcjr12RCTxxVT6xpgA8BegDxACfgWsBh4ketLT\nldbamZ0TUY6UHTvK+fVN99PUHGbE8L4MGTKAJ2evxgl2B+CF+VsZ2P8tJpxyssdJRaSzxDrSvxTY\nY629zBjTDVje8e96a+0bxpi7jDHnWmuf7bSkclhc1+Vvf/sH27ftof+AUs47bwbf+s7v2VNTjBsJ\n89biFwgGm8gtnfDJ1/iCOaxbv0mlL9KFxFr6TwBPdlz2Ae3AKGvtGx3LXgKmACp9j7W0tPDkE7OY\nP38RtXU9CfjTWbDgDZ548nnqmgvJyII9lWso6H4crS311O7dTE5+XwCccBWjR53k8SMQkc4U04Zc\na22jtbbBGJNFtPx/DOw/8VsH5HRCPjkMbW1tfPfaX/L8M+Vs2dxOwJ9OQ+Me2sLthN2hNDXWAOD4\nfPj9IdLS8/C5fvZsf5fSbtVc/oUTGDniOI8fhYh0ppg35BpjyoCngT9aax83xvxuv5uzgOpDWU9h\nYVasEeJCPOd/7tmX2FWeTSAQwHWjx8ZpbKykoPgYAFKCmewqX4brNuC6ERzHR3ZuGYW5KTz+0K+8\njH7I4vn5PxSJnD+Rs0Pi549VrBtyi4GXgZnW2tc6Fi8zxkyw1i4EpgHzD2VdlZV1sUSIC4WFWXGb\nv62tjU2bt+ESLfu0lG7sqVxDa1sDkUgYn89PTk4vUkIhfnHjhfz1iTlsL68jJyvI16/6Qtw+rv3F\n8/N/KBI5fyJnh66RP1axjvR/BOQCPzHG/BRwgW8BdxhjgsAa4KmYU8lhqaio4KfX387ucthXt4Gi\nwhPIzCimKKuNK6+6hEcenU35zhA+Xzunndqb4cOHM3z4cK9ji8hREFPpW2u/DXz7ADdNOqw00inu\nu/tx6qqKSE91SEspZF/9Uq646gJmnHU1aWlpjBkzhnXrLH36dCctrZvXcUXkKNKHs7qQD5av4P4/\nP826dZsf83gqAAALnklEQVTJyx4BEJ2nz+zBWWdPIzU1FQC/38+QIcck/FtcEfnvqfS7gO3btvPc\nsy/z8guLyQwdgxPJpqGxkoz0QiKRdvoMyPyk8EUkuan0E9yO7Tv48XfvoLk6m9aWNAhBXnZf9tVt\noSVSyamTT+TKq/4/r2OKSJxQ6Se452a9Qri+mIDfpaWl9pPlWenFTJ5RypVX6UBpIvIpHWUzga1d\na1m9chURN4zjOORm9GTn7mUEM8sZe2ou/3vlpV5HFJE4o5F+gpr93Ms8ds8CIi15VNS8R3Husfh9\nPiacfgw/++X3dWRMETkglX6CWbvW8viDz7N0yRpyQkPw+aE0dyRtqev57vVfY8yYE1X4IvIfqfTj\n3PJly9m4/iPKevXkvcXLmD97OWmRvrTUA3nR+/h8AQryShk7doynWUUk/qn049gDdz/CK39bSVNj\nC61OAz6CZKSVQAD8TpD6pt1kphUR9tUyYcoIr+OKSAJQ6ccJ13V59P6/sWX9Top6duMrV13Kq7Pe\nIxgppqp9HcW5g2lqraaxpYpQoCcFWf3Y17CVvKG7+eIl5zPuZB0CWUQ+m0o/Ttx56z0s+vtmAqSy\nOrKPyoo/4EY6bnRdANJCuTTU72NnzUoKCvIZO6UPP77xu5rDF5FDptKPE+uXbydA9FOzrhtm0Wvv\nUtK7gKbtTWSlFLG7bi2FmYPIyShgzNTufOO7VxMI6McnIv8dtYbHPlz2IY/d/TQb1m0kD0NbuJnK\n5k10Dw6jZVOYttyNTJ46kcHHzqBi5y569CxlwqRTvI4tIglKpe+hurpabv3x/YR3ZZPans+OtlWE\nnTbKskfiOA5+J0Coui89epdw6umTvI4rIl2ASv8oq6qqYtvWrQwcNAi7xtKwwyE1AOmBXNL82bi9\ny3ErXJyOs0+6boRQKOhxahHpKlT6R9ELT8/mr79/nrZqH/7ujZjj+9KeWgft0bPghJ12pp11Om+/\n/gE169OIuGF6jvYxZeoUj5OLSFeh0j8KXNelomInT9z1HIG6XFrcaup3tLFpbyvNLW20ZmygsLCI\nY8f24dKvXsKFl3yeObNfIRQKcub0M/D7/V4/BBHpIlT6R5jrutx43S9ZMWcTda21lKTmUh+upjhj\nAAAFKWU0uXu56YEfUVxcAkBqairnXXCOl7FFpItS6R9BNTXV3HPHPdgXd5NON6rDe2kLtxA9pfCn\nHJ+L368fhYgceWqaTua6Lq/NfY1l777PspfXsnPrLrqnDgQHuof6UdmyjcKhGTRX7yRUW0TY18yJ\nZw2moKDA6+gikgRU+p2ooaGBW2+4jRVPb6KqbTfdU/uTH/BT2bqFopQ+OI5D9/6F/PbenxEIBJj/\n8muU9CjhlInjvY4uIklCpd8JWltb+cnMn7LmrQ1E6v3k+grxOdHz06T40slyu7E7sJ7J50/iwss/\nR3FxMQAXfunzXsYWkSSk0u8ED/zxAbbN2UfQTafN1waA3/VT315NZiCXYCCFKRdM4robvuNxUhFJ\ndir9GLwxfyHPPzCb9nALtfX17NpaSY7Tg1Qy2Bv+iDQnk3x/d2q67aT/pL4cM2IIF16qUb2IeK9T\nS98Y4wB3AscBzcAV1tqPOvN7eO2j9Ru569r7oCLATrbS3dcHFz91vn1k+bvR3elLXXE5Z3/pLKad\n/33KevfyOrKIyCc6e6R/HpBirR1njBkD3NKxrMt4a8Ei3J1+HAcCThDHccghn+rIHmryyxl7xolc\nfu1PKe1R6nVUEZF/4+vk9Y0H5gBYaxcDx3fy+j1nhg0hktkKQNht/2R5DvmMnjSCH//+ehW+iMSt\nzi79bKBmv+vtxpjO/h6eOmHsCcz43hkEB0Yo7J1HS699pAyEXtML+M4vr/U6nojIQTmu6372vQ6R\nMeZm4G1r7VMd17daaw82qd1531xEJHnEfLq8zp7TXwScBTxljBkLrPisL6isrOvkCEdPYWGW8ntI\n+b2TyNmha+SPVWeX/jPAFGPMoo7rl3fy+kVE5DB0aulba13gms5cp4iIdJ4utZFVREQOTqUvIpJE\nVPoiIklEpS8ikkRU+iIiSUSlLyKSRFT6IiJJRKUvIpJEVPoiIklEpS8ikkRU+iIiSUSlLyKSRFT6\nIiJJRKUvIpJEVPoiIklEpS8ikkRU+iIiSUSlLyKSRFT6IiJJRKUvIpJEVPoiIklEpS8ikkRU+iIi\nSSQQyxcZY7KBR4FsIAh8x1q72BgzFrgNaAPmWmt/3mlJRUTksMU60v8O8Kq1dhJwOXBnx/K7gC9a\na08BxhhjRhx+RBER6SwxjfSBW4CWjstBoMkYkwWErLWbO5a/DJwOLD+shCIi0mk+s/SNMV8FrgVc\nwOn4/3Jr7VJjTAnwCPBNolM9tft9aR3Qt9MTi4hIzD6z9K21fwH+8q/LjTHHAn8FvmutfbNjpJ+9\n312ygOrOCioiIofPcV33v/4iY8wxwD+AL1hrV+y3/H3gAmAz8AJwg7X23c6JKiIihyvWOf1fAynA\nH4wxDlBtrT0fuIbo6N8HvKLCFxGJLzGN9EVEJDHpw1kiIklEpS8ikkRU+iIiSUSlLyKSRGLdeydm\nXeG4PR17LN0JHAc0A1dYaz/yNtXBGWMCRD9v0QcIAb8CVgMPAhFgpbV2plf5DoUxpgh4D5gMhEms\n7D8EziH6mr8TWEiC5O947TxE9LXTDlxJgjz/xpgxwE3W2lONMf05QGZjzE+BGUS759p42uvwX/KP\nAG4n+jNoAS6z1lYaY64EvkY0/6+stbMPtk4vRvpd4bg95wEp1tpxwI+IHpYi3l0K7LHWTgCmAX8k\nmvt6a+1EwGeMOdfLgAfTUTx/Bho7FiVS9onASR2vl0lALxIoPzAd8FtrTwZ+QXSX7bjPb4z5HnAv\n0d3L4QCZjTEjgQnW2jHAxcCfvEn77w6Q/zZgprX2NOAZ4AfGmGLgG8BJwFTgN8aY4MHW60Xp3wLc\n3XH5s47bE6/GA3MArLWLgeO9jXNIngB+0nHZR3S0MMpa+0bHspeIjqDj1e+JDgzKiR4OJJGynwms\nNMbMAp4j+sHFRMq/Dgh0vMPNITqiTIT8G4Dz97s++l8yTyH6u/wKgLV2G+A3xuQf1ZT/2b/mv2i/\nD8MGiM4ynAi8aa1tt9bWAuuB4Qdb6REtfWPMV40xK4wxH378PzDQWtuy33F7fsiBj9uTcySzHaZs\noGa/6+3GmLjePmKtbbTWNnT8gX0S+DHR8vxY3D7nxpivALuttXP5NPP+z3fcZu9QAIwGPk/0A4yP\nkVj564keR2st0QHb7STAa8da+wzRwc3HDpQ5i3/+Xa4nTh7Lv+a31u4CMMaMA2YCt/LvXfSZ+Y/o\nnH4XPm5PLdGMH/NZayNehTlUxpgy4Gngj9bax40xv9vv5nh+zi8HIsaYKUS3ozwMFO53ezxnB6gC\n1lhr24F1xphmoOd+t8d7/muBOdbaHxtjegCvE90u9LF4z/+x/X9Hs4B9RH+XE6Z7jDEXEZ1Snm6t\nrTLG/Nf5j/rotOO4PU8AX7LWfvy2qg5oMcb07XgLeSbwxkFW47VFROc56dgAveLgd/dex9zfy8D3\nrbUPdSxeZoyZ0HF5GnH6nFtrJ1prT7XWnkr0UN1fBl5KhOwd3iQ634oxphTIAOZ1zPVD/Offy6ej\nyWqig8VlCZT/Y+8f4DXzFnCGMcYxxvQCHGvtXs8SHoQx5lKiI/xJ1totHYuXAOONMSFjTA4wGFh5\nsPUc9b136BrH7XkGmGKMWdRx/XIvwxyiHwG5wE869lZwgW8Bd3Rs+FkDPOVhvv/WdcC9iZDdWjvb\nGHOKMWYJ0SmGa4gelPC+RMhPdAPiX4wxC4luh/shsJTEyf+xf3vNWGtdY8wbwNtEfzbxuheSD/gD\nsAV4xhjjAgustTcaY24nOrBwiG6obj3YunTsHRGRJBLXGx9FRKRzqfRFRJKISl9EJImo9EVEkohK\nX0Qkiaj0RUSSiEpfRCSJqPRFRJLI/w/yvHXRscvyuQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x114d28a90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# the color variation shows the sequence order\n",
    "plt.scatter(a[:,0],a[:,1],c=c_space,cmap=cmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        ,  0.99951139],\n",
       "       [ 0.99951139,  1.        ]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check the the correlation between the sequences (off-diagonal elements), \n",
    "# which will be high for small noise\n",
    "np.corrcoef(a,rowvar=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# now create two-similarly related sequences, but with non-even spacing and larger noise\n",
    "_ = [1,4,5,8,15,45,48,50,55,60,88,89,90,93,99]\n",
    "b = np.array([[i,i+np.random.normal()*5] for i in _])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x114e399b0>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAECCAYAAAASDQdFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF3BJREFUeJzt3Xl8VPW9xvHPTFaME8AQUASL2/26XPeqBRXQioob7tpq\ne91vlVutrVqXarWr1/pCi9blatWuWkqlanFXlKXuoIKlX7e6oICpGAg7Seb+MYMOiCGemeTMzO95\n/+OcZQ6PYfLkx+8sSaTTaUREJAzJuAOIiEj3UemLiAREpS8iEhCVvohIQFT6IiIBUemLiASksrM7\nmtmewFXuvq+Z7QyMBVqBFcA33b3JzM4AzgRWAT9194ldEVpERKLp1EjfzC4AbgVqsquuA0a7+37A\nBOD7ZtYP+DYwGDgI+LmZVRU+soiIRNXZ6Z03gCNzlo9395nZ15XAcmAPYKq7t7r7IuB1YMeCJRUR\nkbx1qvTdfQKZqZzVy/MBzGwIMBq4FqgHFua8bTHQs2BJRUQkb5FP5JrZ8cCNwMHu/hGwiEzxr5YC\nmvOLJyIihdTpE7m5zOwkMidsh7v76mJ/DviJmVUDPYBtgFkdHSedTqcTiUSUCCIiIYtcnF+49M0s\nCfwSeAeYYGZp4Cl3v9LMxgJTs4EucfeVHR0rkUjQ1NQSIXZxaGxMKX+MlD8+pZwdyiN/VJ0ufXd/\nBxiSXWz4nH1+Dfw6choREelSujlLRCQgKn0RkYCo9EVEAqLSFxEJiEpfRCQgKn0RkYCo9EVEAqLS\nFxEJiEpfRCQgKn0RkYCo9EVEAqLSFxEJiEpfRCQgKn0RkYCo9EVEAqLSFxEJiEpfRCQgKn0RkYCo\n9EWk7C1sXsCkx+/ipemT444SO5W+iJS1eXPf5ZlJ3+CwwVezxUbnMXHCFXFHipVKX0TK2ssv3MbX\nD59HTU2SzQfC5htPZP78+XHHio1KX0TKWkVF+xrLtdVtrFq1KqY08VPpi0hZ28KO44FJvQBYtDjN\njNf2YtNNN405VXwq4w4gItKVtthyR6qrb+buRx6gqnojjjz+RBKJRNyxYtPp0jezPYGr3H1fM9sS\nuBNoB2a5++jsPpcDhwCrgPPc/fnCRxYR+WIGDNyaAQPPjTtGUejU9I6ZXQDcCtRkV40BLnH3YUDS\nzEaZ2S7AUHffE/ga8KuuCCwiItF1dk7/DeDInOXd3H1K9vWDwAhgb+ARAHd/D6gws4ZCBRURkfx1\nqvTdfQLQmrMqd0KsBegJpICFOesXZ9eLiEiRiHoiN/caqBTwMbAIqF9rffP6DtTYmIoYoTgof7yU\nPz6lnB1KP39UUUt/upkNdffJwEjgCeBN4H/N7BpgIJBw9wXrO1BTU0vECPFrbEwpf4yUPz6lnB3K\nI39UUUv/fOBWM6sCZgPj3T1tZlOAp8lM/4yOnEpERLpEIp1Ox/nnp0v9p63yx0f541PK2aEs8ke+\n0UB35IqIBESlLyISED2GQURKxvSXJvPB3FlsveVg7D92iTtOSdJIX0RKwsOP3syKqnPYdfgNzPn4\nTJ5+9t64I5Uklb6IFIXFi1tYsWLF527/eOkEvrTFSgC23XEJc+aP665oZUXTOyISq7a2Nn7xh8v5\nd883Sa+oYI/6ERx/0Mmf3TG95gUr6XS4T8rMh0b6IhKrPz/yO1r3ep/G3Tek7949eDHxCO+8+/Zn\n9tsodRT/er2GdDrNqy9tyGabnND9YcuARvoiEqvFqxZS1ePTKqrdJMHcD+fwpc0GrbHfAV89k1dm\n/icvTZ6JbT2ErbbcoZuTlgeVvojEZsGCj6hpraf51RX02j7z5PaVL9eyw1HrvjJnxx2GsOMOQ7oz\nYtlR6YtILJ564Slu83Gs2qySVc+uZNv3etO3vh9f/+pJ1NXVxR2vbKn0RSQW41+dSMUevagAag/Z\nimUvwDnHXBp3rLKnE7kiEovWZPsay6uSbTElCYtKX0RiscMGW9G6YDkArR8uY+f6bQBIp9OsXLky\nzmhlTdM7IhKLs4/+Fv0n3cvbs+ewdZ/NOeSwg5n63ESe+detVG6wHJq35L+PuYYePXrEHbWsqPRF\nJDZH7Dvqk9fLly/nuTm3sN2BmbtyW1e9xl8evp6TjrgwrnhlSdM7IlIUFi5spkfD4k+WK6uStFUs\n7OAdEoVKX0SKQmNjX5a8N4DVv9jpoznt9O+1c8ypyo+md0SkKCSTSf7roGu4f9JY0pVL6b/hbozY\n79i4Y5Udlb6IFI3GPn059aifxB2jrGl6R0QkICp9EZGAqPRFRAKi0hcRCYhKX0QkIJGu3jGzSuA3\nwCCgFTgDaAPuBNqBWe4+ujARRUSkUKKO9A8GKtx9L+DHwM+AMcAl7j4MSJrZqI4OICIi3S9q6b8G\nVJpZAugJrAJ2dfcp2e0PAvsXIJ+IiBRQ1JuzFgObA/8EGoDDgH1ytreQ+WEgIiJFJGrpnwc85O6X\nmtmmwJNAdc72FNDcmQM1NqYiRigOyh8v5Y9PKWeH0s8fVdTSX0BmSgcy5V4JzDCzYe7+FDASeKIz\nB2pqaokYIX6NjSnlj5Hyx6eUs0N55I8qaulfB9xuZpOBKuAi4EXgNjOrAmYD4yOnEhGRLhGp9N19\nCXD8OjYNzyuNiIh0Kd2cJSISEJW+iEhAVPoiIgFR6YuIBESlLyISEJW+iEhAVPoiIgFR6YuIBESl\nLyISEJW+iEhAVPoiIgFR6YuIBESlLyISEJW+iEhAVPoiIgFR6YuIBESlLyISEJW+iEhAVPoiIgFR\n6YuIBCTSL0YXkWiWLl3Kz8bdxPz0YgbW9OS7h5/Ghhum4o4lAdFIX6Qb/fjuG/j71q28te0GPDlo\nJVfefUPckSQwKn2RbjQn3UKiIvNtl0gmeJ/FMSeS0Kj0RbpRn3Qt6XT6k+WGdG2MaSREkef0zewi\n4HCgCrgRmAzcCbQDs9x9dCECipSTiw4/nR/99RbmsYSB1fVccNhpcUeSwEQqfTMbBgx29yFmVgec\nD4wBLnH3KWZ2k5mNcvd7CxlWpNT169uPX515OQCNjSmamlpiTiShiTq9cyAwy8z+CtwH/A3Y1d2n\nZLc/COxfgHwiIlJAUad3+gCbAYcCW5Ap/twfIC1Az/yiiYhIoUUt/Y+A2e7eCrxmZsuBATnbU0Bz\nZw7U2Fja1ygrf7yUPz6lnB1KP39UUUt/KnAOcK2Z9QfqgMfNbJi7PwWMBJ7ozIFKeU6z1OdklT9e\npZy/lLNDeeSPKlLpu/tEM9vHzJ4DEsBZwNvAbWZWBcwGxkdOJSIiXSLyJZvuftE6Vg+PHkVERLqa\nbs4SEQmISl9EJCAqfRGRgKj0RUQCoufpi+TpwSmTeGXOW2zR0I9jRhxCIpGIO5LI51Lpi+ThjvvH\nc0uz09anF3z8D97+4zwuOPH0uGOJfC5N74jk4cn3X8sUPkB9HVMXvBtvIJH1UOmL5KEqveZUTk1a\nUztS3FT6Ink4bfAIer/+AW2LFlP35jy+ucvQuCOJdEhz+iJ5GLzzl/nTl7bg5X++yrbDt6Zfv43j\njiTSIZW+SJ56996I4YP3iTuGSKdoekdEJCAqfRGRgGh6R2Qtf3zgPmZ9MIeN6+o4+9gTqazUt4mU\nD430RXL831/uZuyc13m8rpLfrWzh0pvHxh1JpKBU+iI5nv/gPdKpzG8lSlZXM7Pl45gTiRSWSl8k\nxwaJ5FrLFTElEekamqyUoE1+9mkem/EyNckEZx93HN8+9Cjev+sO3q2uoGH5Ks7ab2TcEUUKSqUv\nwZr2/PNc+fAkVjb0Jd2aZva113HHD37AXedfxvz582ho6EOPHj3ijilSUJrekWA9+fJLrGzoC0Ai\nkeDNyh68/sbrVFVVMWDAQBW+lCWVvgSrrrqKdFvbJ8u1K5fTsFFDjIlEup5KX4L1reNOYLuF80l/\nOJfque9y4o7b0a9fv7hjiXQpzelLsGpra7n50h8wf/486urqqK/vGXckkS6n0pegJZNJNtmkf9wx\nRLpNXqVvZn2BF4D9gTbgTqAdmOXuo/NOJyIiBRV5Tt/MKoGbgaXZVWOAS9x9GJA0s1EFyCciIgWU\nz4nca4CbgA+ABLCru0/JbnuQzOhfRESKSKTSN7OTgQ/d/VEyhb/2sVoAnRUTESkyUef0TwHazWwE\nsBPwW6AxZ3sKaO7MgRobUxEjFAflj5fyx6eUs0Pp548qUuln5+0BMLMngG8BvzCzoe4+GRgJPNGZ\nYzU1tUSJUBQaG1PKHyPlj08pZ4fyyB9VIS/ZPB+41cyqgNnA+AIeW0RECiDv0nf3/XIWh+d7PBER\n6Tp6DIOISEB0R66Utba2Nq697XbebGomVV3B904+iX59+8YdSyQ2GulLWRt7x2+Y+P5KXqc3L65I\ncdkNt8QdSSRWKn0pa/9qWkhFdQ2QeWb++4uW097eHnMqkfio9KWsNWxQTTqn5HvXVpFM6mMv4dKn\nX8ra+aefzE4VzaSa32Pg8rlc+I3j4o4kEiudyJWyVldXx5hLL4w7hkjR0EhfRCQgKn0RkYCo9EVE\nAqLSFxEJiEpfRCQgKn0RkYCo9EVEAqLSFxEJiEpfRCQgKn0RkYCo9EVEAqLSFxEJiEpfRCQgKn0R\nkYCo9EVEAqLn6ZeBJ6dM47FpM0gm03zj6EPYesst4o4kIkVKI/0S9+JLL3PtXZN46cMaps+r5bIx\nd9Lc/HHcsUSkSEUa6ZtZJXA7MAioBn4K/AO4E2gHZrn76MJElI78/bkZtNX2+2S5OdGHZ557gYMO\nGBFjKhEpVlFH+icB/3b3ocBI4AZgDHCJuw8DkmY2qkAZpQN9etfTtmrZJ8vJlYsYtNnAGBOJSDGL\nWvrjgMtyjtEK7OruU7LrHgT2zzObdMLxRx/B7v3bSLS8S+Xitzlqn63YZptt4o4lIkUq0vSOuy8F\nMLMU8GfgUuCanF1agJ55p5P1SiaT/Oji81i2bBmVlZVUVVXFHUlEiljkq3fMbCBwD3CDu99tZlfn\nbE4BzZ05TmNjKmqEolA8+aPlKJ780Sh/fEo5O5R+/qiinsjtBzwMjHb3SdnVM8xsqLtPJjPP/0Rn\njtXU1BIlQlFobEwpf4yUPz6lnB3KI39UUUf6FwO9gMvM7HIgDZwLXG9mVcBsYHzkVCIi0iWizul/\nB/jOOjYNzyuNiIh0Kd2cJSISEJW+iEhAVPoiIgFR6YuIBESlLyISEJV+EZg+fQaPPT6JZcuWrX9n\nEZE86Hn6Mfv5VdfzzPRmEoke3DXuUcZcfSE9e/aKO5aIlCmN9LvQu++8w92//xMzX5m5zu1vvPkG\nT7/4EdU1DVRVb8DCJf357e/v6eaUIhISjfS7yLSnpnHzFX+Cj+u4t3Yao87emxO+edwa+yxdspR0\n+tO/gkQiSVtburujikhANNLvIvf94VESzSkSiSRVK1I8Om7aZ/bZfvvtGbRpK+3trQBUJuYy8sCh\n3R1VRAKikX4XSbevOWJvX8cIvqKigmuuvpjf/2E8y1es4qADTmFL/X5bEelCKv0uMvywr/DH2Q+R\nXLohrRXL2OvAHde5X01NDaedemI3pxORUKn0u8jBo0bSr39fZjz7MoO22oz9D9IvEhOR+Kn087Bg\nwQJ++D8/YeHcFjbdZmPO/eF3qK6u/mT7brvvxm677xZjQhGRNelEbh4uOvUyXhv3Lh9ObebFW1/l\nuiuuizuSiEiHVPp5+GD2hyQSCQCSiQo+8LkxJxIR6ZhKPw8bDaj/5HU6nabXJvUd7C0iEj+Vfh4u\nuf5CGvaqo2LzNgYc1MC5Pz437kgiIh3Sidw8bLOdcf091693vyVLljBh/P2QgKOPHUWPHj26IZ2I\nyGep9LvYkiVL+O7ZV7J4Xh/SwLQnr+TaG6+gtrY27mgiEiBN73SxCePvp2VeHxKJJMlEkoVzenPv\nhIlxxxKRQKn0u1gimQA+fQRDmjTJZCK+QCISNJV+Fzv62FH0HtBMe3srbe2r6DOohVFHHhp3LBEJ\nlOb0u1htbS1jbryC++97gGQiwaGHH7zGXbsiIt2poKVvZgngRmAnYDlwuru/Vcg/oxTV1NRwzLFH\nxh1DRKTg0ztHADXuPgS4GBhT4OOLiEgeCl36ewMPAbj7s8CXC3x8ERHJQ6FLvx5YmLPcamY6WSwi\nUiQKfSJ3EZDKWU66e3tHb2hsTHW0uegpf7yUPz6lnB1KP39UhS79acChwHgz+wowc31vaGpqKXCE\n7tPYmFL+GCl/fEo5O5RH/qgKXfoTgBFmtvq3gJ9S4OOLiEgeClr67p4GzirkMUVEpHB0klVEJCAq\nfRGRgKj0RUQCotIXEQmISl9EJCAqfRGRgKj0RUQCotIXEQmISl9EJCAqfRGRgKj0RUQCotIXEQmI\nSl9EJCAqfRGRgKj0RUQCotIXEQmISl9EJCAqfRGRgKj0RUQCotIXEQmISl9EJCAqfRGRgKj0RUQC\nUhnlTWZWD/weqAeqgO+6+7Nm9hXgOmAV8Ki7/6hgSUVEJG9RR/rfBR5z9+HAKcCN2fU3ASe4+z7A\nnma2c/4RRUSkUCKN9IExwIrs6ypgmZmlgGp3fzu7/mHgq8BLeSUUEZGCWW/pm9mpwHlAGkhk/3uK\nu79oZhsDvwPOITPVsyjnrS3A5gVPLCIika239N39duD2tdeb2Q7AH4HvufvU7Ei/PmeXFNBcqKAi\nIpK/RDqd/sJvMrPtgL8Ax7n7zJz104GjgbeBvwFXuPvzhYkqIiL5ijqn/zOgBvilmSWAZnc/EjiL\nzOg/CTyiwhcRKS6RRvoiIlKadHOWiEhAVPoiIgFR6YuIBESlLyISkKhX70RWDs/tyV6xdCOwE7Ac\nON3d34o3VcfMrJLM/RaDgGrgp8A/gDuBdmCWu4+OK19nmFlf4AVgf6CN0sp+EXA4mc/8jcBkSiR/\n9rPzGzKfnVbgDErk629mewJXufu+ZrYl68hsZpcDh5DpnvOK6arDtfLvDIwl83ewAvimuzeZ2RnA\nmWTy/9TdJ3Z0zDhG+uXw3J4jgBp3HwJcTOaxFMXuJODf7j4UGAncQCb3Je4+DEia2ag4A3YkWzw3\nA0uzq0op+zBgcPbzMhzYjBLKDxwMVLj7XsCPyVyyXfT5zewC4FYyl5fDOjKb2S7AUHffE/ga8Kt4\n0n7WOvJfB4x29/2ACcD3zawf8G1gMHAQ8HMzq+rouHGU/hjgluzr9T23p1jtDTwE4O7PAl+ON06n\njAMuy75Okhkt7OruU7LrHiQzgi5W15AZGHxA5nEgpZT9QGCWmf0VuI/MjYullP81oDL7L9yeZEaU\npZD/DeDInOXd1so8gsz38iMA7v4eUGFmDd2a8vOtnf/4nJthK8nMMuwBTHX3VndfBLwO7NjRQbu0\n9M3sVDObaWavrP4vsLW7r8h5bs9FrPu5PT27Mlue6oGFOcutZlbU50fcfam7L8n+gP0zcCmZ8lyt\naL/mZnYy8KG7P8qnmXO/3kWbPasPsBtwDJkbGP9AaeVfTOY5Wv8kM2AbSwl8dtx9ApnBzWrrypxi\nze/lxRTJ/8va+d19PoCZDQFGA9fy2S5ab/4undMv4+f2LCKTcbWku7fHFaazzGwgcA9wg7vfbWZX\n52wu5q/5KUC7mY0gcx7lt0BjzvZizg7wETDb3VuB18xsOTAgZ3ux5z8PeMjdLzWzTYEnyZwXWq3Y\n86+W+z2aAj4m871cMt1jZseTmVI+2N0/MrMvnL/bR6fZ5/aMA77u7qv/WdUCrDCzzbP/hDwQmNLB\nYeI2jcw8J9kT0DM73j1+2bm/h4EL3f032dUzzGxo9vVIivRr7u7D3H1fd9+XzKO6vwE8WArZs6aS\nmW/FzPoDdcDj2bl+KP78C/h0NNlMZrA4o4TyrzZ9HZ+ZvwMHmFnCzDYDEu6+ILaEHTCzk8iM8Ie7\n+zvZ1c8Be5tZtZn1BLYBZnV0nG6/eofyeG7PBGCEmU3LLp8SZ5hOuhjoBVyWvVohDZwLXJ898TMb\nGB9jvi/qfODWUsju7hPNbB8ze47MFMNZZB5KeFsp5CdzAvF2M5tM5jzcRcCLlE7+1T7zmXH3tJlN\nAZ4m83dTrFchJYFfAu8AE8wsDTzl7lea2VgyA4sEmRPVKzs6lp69IyISkKI++SgiIoWl0hcRCYhK\nX0QkICp9EZGAqPRFRAKi0hcRCYhKX0QkICp9EZGA/D/dk8ccOjj9FgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x114d560f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(b[:,0],b[:,1],c=np.linspace(0,99,len(b)),cmap=cmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        ,  0.99154274],\n",
       "       [ 0.99154274,  1.        ]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# check the correlation\n",
    "np.corrcoef(b,rowvar=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# now create randomly-ordered seqences with larger noise\n",
    "\n",
    "\n",
    "_ = np.array([random.random()*100 for _ in range(100)])\n",
    "c = np.array([[i,i+np.random.normal()*10] for i in _])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x114f3d278>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAECCAYAAAASDQdFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XdgFNXax/HvbEmvpACB0AROQECaIB2VIlgoYkcvYL9c\nr3rt2LuvBbtiw66oFBtFVDqCgEiHQ5ESaoB0Und33j82hAQCJJslyWafzz9kdmfO/HbZPHtyZuaM\nYZomQggh/IOlugMIIYSoOlL0hRDCj0jRF0IIPyJFXwgh/IgUfSGE8CNS9IUQwo/YyruiUqor8ILW\n+nylVHvgDcAB5AM3aK0PKqVuBm4BCoFntdbTz0RoIYQQnilXT18pdR/wARBY9NBrwFit9QXANOAB\npVRd4A6gG3AR8LxSyu79yEIIITxV3uGdrcCwEstXaa3XFv1sA/KALsAirbVDa50JbAHaeS2pEEKI\nSitX0ddaT8M9lHN0+QCAUqo7MBZ4FYgAMkpslg1Eei2pEEKISvP4QK5S6irgHWCw1vowkIm78B8V\nDqRXLp4QQghvKveB3JKUUiNxH7Dtq7U+WtiXAc8opQKAYCAJWHeqdkzTNA3D8CSCEEL4M48LZ4WL\nvlLKArwO7ASmKaVMYL7W+kml1BvAoqJA47TWBadqyzAMDh7M8iB2zRAXFy75q5Hkrz6+nB1qR35P\nlbvoa613At2LFmNOss5HwEcepxFCCHFGycVZQgjhR6ToCyGEH5GiL4QQfkSKvhBC+BEp+kII4Uek\n6AshhB+Roi+EEH5Eir4QQvgRKfpCCOFHpOgLIYQfkaIvhBB+RIq+EEL4ESn6QgjhR6ToCyGEH5Gi\nL4QQfkSKvhBC+BEp+kII4Uek6AshhB+Roi+EEH5Eir4QQvgRKfpCCOFHpOgLIYQfkaIvhBB+RIq+\nEEL4EVt5V1RKdQVe0Fqfr5Q6C/gEcAHrtNZji9Z5DLgYKATu1lov935kIYQQnipXT18pdR/wARBY\n9NB4YJzWug9gUUoNUUp1AHprrbsC1wBvn4nAQgghPFfe4Z2twLASy5201guLfp4J9Ad6ArMBtNbJ\ngFUpFeOtoEIIISqvXEVfaz0NcJR4yCjxcxYQCYQDGSUezy56XAghRA1R7jH947hK/BwOpAGZQMRx\nj6efrqG4uHAPI9QMkr96Sf7q42vZUw4e4OOfX8NpzSWpbheGD7y2uiNVC0+L/kqlVG+t9QJgEDAH\n2Ab8n1LqZSARMLTWqadr6ODBLA8jVL+4uHDJX40kf/WpjuyzFs/mn5RdtGvcmp4du1doW6fTyTOf\n30HC4CwMw2BZ8gZypzkZ0HPoGUp7ZlXmC9fTon8v8IFSyg5sBCZrrU2l1EJgCe7hn7EepxJCiBLe\n+/4jfo3YhK1FCPOSN7H/twOM6Dfs9BsWSUk5gLXJQQwjGIDIRBvb/1gN+GbRr4xyF32t9U6ge9HP\nW4C+ZazzFPCUt8IJIXxbenoai/5aTtOGiZytWnnczrKMjdiahwJgSQxl4aqVjKB00c/Ly+PneTOw\nWmxccv4g7HZ78XNRUdE4DgUWLzsLXQSYvjU85S2e9vSFEOKUNv+zlbunfsK+xFjsm5fyr7Vn8e8R\n13nUlsVllF42Sy/n5uZy58TH2X9uCKbTxW/vL+G1W54qLvzBwcH0bvAv5v32BUZoAVE5zfjX0H97\n9sJ8nBR9IcQZ8dHcGRxo0RAL4EwI4ju9npsLC4sL8ZwlM9mSspZwWzRXXTQGq9V60rYuadaHLzf9\nitk0DNvmbIaffSUfTv6GJTt3EWBAYpiF/V1DsdjcbezqCD/NmcHwgUOK27ig2yX07TqYgoICEhPj\nfPZ4SmVJ0RdCnBGO45adVgOHw4Hdbmf6vMksyP4Me95BHIUm/3y4kUduHX/Stob0uYS2O1qzbus6\nOvXtyJ/r1vNx8m7MqGgA1v41j5hOLYvXN6wWnKaDZUuW8/FbU8jNKqRxUiwPP3MPQUFBZ+Ll+gyZ\ne0cIcUZc1vZcQvccBMDMyaNHaCzBwe4DqSt3zCWxYBs3XZbNrcOzwZxFauqpT/Zr1qQZl/W7jAb1\nG7I+eRdmxLHLgHIbtCRo7n5Ml4nL4aTusiMM7nUR7770NUeSI3Clx7Dtj0ImvD7xzL1gHyE9fSGE\nV5imCYBhuMfb+3TpxvjQUOavXUl8eF2uu/rYmTJZ+5O5/ipL8fpXX2th2dIZXDRg5Gn343K5MFOX\n04pk9uXVJa1eD2IdDl4d/TCzV8zDarEw4sah5OXlciTVxdF+vcViI/VgtndftA+Soi+EqLRZcz9n\n5a5PyXPksG1dDP26XM/oa66k09nt6HR2uxPWb9+sD7k52wgOcRf+9FSTiIjYUuuYpklKygFCQ0MJ\nCwsnNfUwe3Zv4++1U7j6hoWEhRukp+3g/944wHWXjiMhoQGjLjt2oDgoKIg6DQLISXYvF7pyadqy\n8Zl7E3yEFH0hRKVs37GVZOMjzhtqAYJo3zeD91+cQZOGCZzfq2eZ21x1+V189MUGWnRYgqPQwv6t\nFzHyyoHFz+fl5XHH6y+x1nARXFjIoCiTbmouSc0PE2YzObjPJCzcTlS0waU9ghjQo9cJ+zAMgwef\nHcuHr39JbnYhzds05oYb/fMq3JKk6AshKmVnsqZ+c7N4OSrWRnCYkw1b/jlp0bdardx8w4ds3boJ\na5idfleeVTwsBPDu5Emsrh+HYbWSA1gz3qVfn3zAxjXD4eufj9C0pfssIEdh6EmzNWnSmGdeHeeV\n11lbyIFcIUSltG3dla1/hhQv61UOco5E0ablWafczjAMWrRoRbNmzUsVfIAsRyFGiVM4Q0LNUs9n\nZwSxaa3B3BmN6Nj2Hi+8Cv8hPX0hRKXExMQysO1zTP36JfanHeLw3iZcd34f+vTs4XGbfVu35bdF\nc8ivG4dpmqzaGMvh1L3E1IEdyVbq1budpIbD6NUhvtSVt3+smMnqnR+DNZ8wZydGDn/0hC8UfydF\nXwhRaa1aduDhll95rb3eXc7jCaeT39avJgiDsfdPYsXqH8nP20FUdHsGXFR6CgbTNHn5xTdwNvia\nHpcFAJCdMZsZvzfi4n6jvJarNpCiL4SokS7o1oMLuh37a6Hv+aNPuu4nE79i3u//MOa1Y8NAYZEW\n9uQln9GMvkjG9IUQ1Wb+4sVMmjaNQ4cPe9zGO298xKQvZhMa2ICVvx0b6jm426RedFtvxKxVpKcv\nhKgWz78zgdl7siE4nMlLXufl/9xIk0YVO49+5vTZLJixB1d+AEaYBf17D9L2riI4PIfenUbS9yL/\nmzr5dKSnL4SococPH+bXrfsxwqIwrFZSoxvx+U8zK9zOjn/2YLOEERfRkoMH15N+OJ3DW1pwzcDX\nGXLRrWcgue+Toi9ELWCaJgUFBdUdo9wcDgeu486qcZ1k3VNp1z4Jp5GOxWIloU5b6sXV4bW3HqZz\n547eCVoLyfCOED5u6Yq/ePPzH8jIc9IwKpBn7h9LbExMdcc6pfj4eLrGh/Jnbi6WwGDC0vcwfOgV\nFW6nR69uHNh/kEVz12C1GQy98jLq169/BhLXHsbRSZKqienLc1r78j1OQfJXN2/lv/5/T3DY1gBw\n9/g7x+Xy9P13VrrdU/FGdpfLxeSffyY1M4sBPXvQrEkT74Qrh1rw2fH44gPp6Qvhw0zTJDPXAUV3\n/jMMg6w8Z/WGKieLxcKVl11W3TH8jozpC+HDDMOgSWwIpss9Iu7Mz0E1jKvmVKImk6IvhI979oH/\n0qOBk7MjjzC8Qwy3jT79nPRnimmaTJ00jQnjP2DF0hXVlkOcnAzvCOHjwsLCeOTumnGT7xcfH8/q\nacnYCGTBN6sY/Vga/S/uX92xRAnS0xdCeIXL5WLN/C2kFx7igDOZ1Kx03hv/aXXHEseRnr4QwisM\nwyDPkYk9MJLYoHoA5Gcc4efvp3PJ0IvZvXcPn8+ciWmaXN2vf5WerSOO8ajoK6VswKdAE9w3vb8Z\ncAKf4L7GYp3Weqx3IgohfIFhGJzVsRF7/zw2/02gNZQ9O/eTmpbKHRMmcDChMWCw9OOJTLjtdhLk\nnPoq5+nwzmDAqrXuATwNPAeMB8ZprfsAFqXUEC9lFEL4iP+NuwtnWFrxsisok07ntWf6nDmk1Ess\nfjy1fiNmLJhXDQmFp0V/M2BTShlAJFAIdNRaLyx6fibQzwv5hBA+pF79etzzf6NJ7GojoZOVkfcN\npHPXTsRGRUFebvF6rrw8IkJDTtGSOFM8HdPPBpoCm4AY4FKg5J2Js3B/GQghqtgXn3zDwl9XYxgG\n/S/rwuVXlv+P7vz8fD754Xly7bux5Edw1fn3UrduQoX23/ncTnQ+txMAGzav5tOf3iTUFkkvw8mC\nlH1gWDgv2M6IwZdWqF3hHZ4W/buBWVrrh5VSDYB5QECJ58OB9PI0FBcX7mGEmkHyVy/JX9q83xcy\n4+t1WJ1RAEz7eBldzmtLx07nlGv71z55kbBefxEVYAEO8s2c53j632WfgXO67EtXLmLKrhep08Ug\nNcdBQmZzvus1hr1799K7Vy/CwsIq9Nq8zdc/O57ytOin4h7SAXdxtwF/K6X6aK3nA4OAOeVpyMfn\nv5D81Ujyn2jpkjVYnSWKWUEEC+YvJ7FRs3Jtn1a4i7iAY6O+Ryx7y8xYnuy/rPiBOt3dU8QEhNjY\nHbqWe0a9iCMtmHfqT+HOJ0bRvmP5voy8rTZ8djzl6Zj+a0AnpdQC4DfgQWAs8KRSajFgByZ7nEoI\n4ZF27VvjsmUUL5uB6XQ+t0O5t7cXxuJyHpuE0V4Q63kYV+nykp2ajyUzjhB7HZyH6vDFhGmety08\n5lFPX2t9BLiqjKf6ViqNEKJSzu3amRE37WPerBUYhsGgoRegklqWe/uRg+9j4o9PkB+4F0tBJCN6\n3uNxliHdb2DCb48SeW4uR/a7ODA7iCDLsVHgAh+ZGK62kYuzhKhlho24lGEjPDtIGhoayh3XvOSV\nHA3qN+KBIW+xYs0S6tVvwIzWc1m+bR92IwgHubTp0sIr+xEVI0VfCFEpDocDm63sUhIWFk7f7gMA\nUA+15qv637J3ZwrNktpy+dVy/9rqIEVfCOGRHcm7eOzTj9mTX0jdABuPXn0trVqcfCjJMAyuG13W\nqLCoSlL0hRCnlZuby7IVK4iNqcPZrc8G4JVvv2FbXfcdu3YCr0ydwocPPHTCtgcOpPDsax9wIC2f\n2MhA7h97PY0bNarK+KIEKfpCiFNKTU3lrmdeZQ8xWBw5DEz6g/97/H+kOxyl1ktzFJa5/cvvfMr2\n7FiMAIPkXBg/4Utef+7EL4eT+WzqVBbpbQQYMGrQADqfUz2nedYWMrWyEOKUJk6awoGgRtiDw7CG\nx/PL+r3s3r2bFpERuAoLADCdTpqHl33ueHp2IYZx7JauqdkF5d73rLlz+WTNP2y2R7POFs3Tk6aS\nmZlx+g3FSUnRF6IWWLN2HtO+v5Ifpw9hxsxXvNp2oWmWKtpOw072kSOMu+lWhgcH0jEni4utJk/d\nWvaNXBLiQjFd7tMzTdNFQkz559xZve0fXGFRxcuHAiNZt3Gjh69EgAzvCOHzMjLS2bX7YQYMOgjA\ngf2ahQsboJIuZOI3Uyk0TQZ0P5eunTp51P7FfXuw7P0p5IYl4HI6aBXppGWLFhw+fIT7Rt9Y5jaZ\nmRl88tUUHE6TKy/tj3X67xxIyyQ2Moj7/nNbuffdKC4G176dWIJDAQjPy6RFs7M8eh3CTYq+ED5u\nxw6NStoPWAGoW8/F2tVr+N//rScluBGGYbDsyxk8ZQ+gQ7u2FW6/3dln8+ztVmYvWkpIQDCjrh6D\nxXLyQYK8vDzuevhlDjkbYBgWlq6exAsPjaFZ0yYV3vfVlw1h5773WZ68B7sB1/fvTVyc3Pi9MqTo\nC+GDnE4netMmgkOCady4JYuX1iUh4RAAB1MsHDoczF5CsRcNyxSE1Wfu0hUeFX2A1klJtE5KKte6\ni/5Ywv68GOxFc/jk2xsy89f5jL2lSYX3axgGD952a4W3EycnRV8IH1NQUMADtzxM8h9pYHPR5fJW\nDBjxNL/MeAdbQD4Gvenb+3ImrfoaiADA5SgkNCi4SvJFR0WBIxcC3GP3LpeTwAApNTWF/E8I4WO+\n/ngSKYvzCTbCwQErvttCv8suYPiwKaXWu+ycRKb//Q+Fho3WUQajry57/N3bOnXsQPfWi1i8LgXT\nsNEsJo/rrxlTJfsWpydFXwgfk3ckD4thLV62OGwcPpx6wnp3jPkXI9NSycnJISGhQakzcM60R+67\ng01ak5OTQ7u2bU86TYOoevI/IYSPufCSC1g87S/MAyGYpkl4a4MevXuUuW50dB2io+tUcUK3JKWq\nZb/i1KToC+Fjmrdszrj37mbWlF+x2Ayuv+1aQkLkfrOifKToC+GDklonkdS6fGfTCFGSXJErhA9w\nHDfPjRCekp6+qHUOp6Xy5KQPOGDm0cAawpPX3kp4eER1x/LI1z9/w5RN32GPNojPj2HcFU9SN65u\ndccSPkx6+qLWefzr91jSLITtzWNY2DiQJ756r7ojeWSt3shnGybRbEQCiRfWJ2CQnXenv17dsYSP\nk56+qHX2mnkYhnvGR8NiYY8rt9Jtulwu3n55AtvW7iEkPICb77mBpmc1rVAbaempfDBzAjnWXBoG\nJPLAjXefcv3fVy3DXsdevGwYBvtzUzzKL8RR0tMXtU49IwjTNAEwTZP6RuWvRJ347qcs+Wo7qetg\n95IC/u/BN4r3UV4vTXuBfV1SyOyczeqmq3n7m3dOuX5sSDhpB0xcThcABdkF1LckePwahADp6Ysa\nZtPW9cxdOwuLaeXK82/w6Bzzx6+4iae++5D9RWP6j119S5nrORwOPps6hay8PAZ2605Si5PfqHv3\nlv3YjIDi5bTduWRlZRIREVmuTC6Xi0MBhwgz3OvbwwPYsXPXKbe59uKhLP9nI0s/XklIuEEjRyzj\nH3qiXPsT4mSk6IsaY8v2TXy05iUiuwRgmibPT32IJ699jeDgivXU68bF8fa/T31nJtM0uevFF/k7\nOAqLPYBfvviKF64aQbuiWwEer079SLaY6VgN969MWFxAhQ4OWywWgguPnUtvukxCXac+t95isfD6\nnQ+TlpaK3R5AWFhYufcnxMnI8I6oMRatm0NkF3dv2jAM7F3zWL5qyRnZ165dO1mZ78Jid+8vK64B\nPyxafNL1b7v7Jlr2jySgQQ6RrRz857HRGIbBr0t+5ZFvHueRrx/nj79PnfVfXUbjWlBI9pJMghcE\n8L8r7ipX1ujoOlLwhdd43NNXSj0IXAbYgXeABcAngAtYp7Ue642Awn8EGEE4CpzYAtzzyhSkOqkT\nF3tG9hUUFITN6eToqLxpmlhPMTVNQEAAj780rtRjf29YyaTU77F3dBfkj1Z9QcP4BBo1aFxmGx1b\nd6Jj606YRXeiiooK5+DBLG+8HCHKzaOevlKqD9BNa90d6As0AsYD47TWfQCLUmqI11IKv3DFgJG4\n5kSQtj2LwxuyaJrSgTat2p2RfdWtW4/Bjetjph3CmZdL/QO7uHHI0Aq18feWVdhbHuuB29qFsmT1\n0uLlXxbNZ8wbLzDqjReY+tus4scrM/HZ6o2reeDzJ7jn60f5+pdvPG5H+C9Pe/oDgXVKqe+BcOB+\n4Cat9cKi52cC/YEfKh9R+IuAgACeGPMKW7dtJjgohMTERh61M/7jiczfthOrASO6dOLaSy8rc70H\nbroZPnyLXakpDBt0IXXj4yu0n0ZxiSw6sIqAuu5jDgU7jpDUxD01gt62hRf+mk92A/eFVNs2r6RR\nfD06t2vv0WsCyM7O4tU/P8TV1X0w+Kd9fxK9OIqLegz0uE3hfzwd048FOgEjgNuBL49rKwso32kN\nQpRgsVho2SLJ44L/06+z+X5/BmnxiRyKS+TDFWvYtGVzmes++clbTIvPZWXneJ7eNp8f588u1z7y\n8vLIycmhX/f+dD3YBtfSbMylRxhAT85pfQ4AS1avJCvh2JdIXnwMyzet9+g1HbVp6yZymh77NbPV\nD2HjgS2ValP4H097+oeBjVprB7BZKZUHNCzxfDiQXp6G4uLCPYxQM0j+6nV8/pTMVAg99lhBZAy7\n9++iV/fSNwV3uVwszdmL0ch93ntBgzrM372RG+MuB2DZ8t/ZvH0RkWGNuGTQqOIhmacnvsO0fWtw\nGdAvoikvjb2/zOGa8887l4++n0RefAwA9rQMzuvV5YS8FXn/z+3QjsCvnVD0XeLILKBF3cRq+z+s\nbZ8df+Fp0V8E/Bd4VSmVAIQCvyul+mit5wODgDnlaciXD2TFxfn2gbjamL914+bY1/9GYZT7AHBE\nWgqtWww9YT3TNLE6Sl9c5cpzcPBgFouWTCPL/gJNuxSQkQZvTFjFtSOeYfGKJXxVuA1a1wfg58w0\n1DffMPTCi0/I1jihGbc1ac2UDaswgUHNkuiQ1KFUjrLy7927h5TUQ7RqkURgYOBxrQYysvFlTFk6\nk0Kri3b2plx89WXV8n9YGz87vqQyX1geFX2t9XSlVC+l1DLAwD3EswP4UCllBzYCkz1OJYSHunXu\nzH8OH+aX1WuxGXDtkMEk1Kt/wnqGYXCN6sp7//xFXnwEdXdlMGbAdQDsOTSTcy4oACAyGhyBizBN\nk+37dmPGhHG0X29EBLM/5cQ7Vh113eChXDe4/AeHP/zpS77L/IvCqAAaznfxytX3ExdT+uylgd0G\nMLDbgHK3KcTxPD5lU2v9YBkP9/U8ihDeMWzgQIYNPP3BzZEXDaPnzg5s2bWdzr3PKb761+m0l1rP\nWRiAYRj069qbb6e+Snob9/hK6OZDnN/3Uq9kzs7OYmrKcoxz6hEAHIg3+eCXrxh37X+90r4QR8kV\nucKvNWnchCaNm5R6rNPZN/Hn75tJ6nSAPf8E0yjmegDqxdflhQtH8cXSWZgGDO18OarZyaduqIgj\nR3IoDLFy9M63hmFQYKnY3D5ClIcUfSGOo1p0ID72GzZsWkanxopGjY7Nptm6RRLPtTj9HatM02Tq\n79NJTj9Ap2Zn06vjeadcPz4+HpUazOZCJxa7Fdu2NC5o0bvSr0WI4xkVnSnQy0xfP5gi+atPTc7/\nwpfvMLPOYYzIEGy70xkb141h5w8utc7x+fPz8/no569Iy88med0+rEY0cZGB3P+fm2rcPXBr8ntf\nHrUgv8dX+MncO0KcAX9k7sCIdBdqR8Moftu5utTzDoeD5ya+z30T3+DdyV/icrkIDAzk35eP5sj2\nQpJzG7MzL4rl+wJ55tUJ1fESRC0lRV/4vfz8/ArPjX86NrN0R8x6XPNPfPIOEwpSmBdv5yPXAV78\n8qPi5/am5mBY3L+ahsXC3tQcr2YT/k3G9IXf2p+yn+d/eIWDwZmE5QVya/dRdGrd0SttX9miOx9u\nW0J+/Qiit2dyQ7erSj2//kgqlnruKRqMkGDWHDx2R6zoUDuHjhxbt05Y6bOJhKgMKfrCJ+zelcyU\nz78H4MrRI6ifcOK59xU1YfZHZPa0EWTE4AAm/vGF14r+1f2H0HVHOzbt2ErnYR2Iiy19vn24YS21\nHMax5Xtvu57n3/qEAxl51I0M4p5bR3klkxAgRV/4gJQDKTx683M4doVgmiarFj7Fy18+U+nL6HOs\neaWmUDhiyz/l+nl5eXz241QKnU5GXDjwtBO0NW3SlKZNyr6P7n8vHMLzsyez12aSkGdy17Dri59r\nkJDAW8+NK3M7ISpLir6o8Wb+MJPCncEYhvv89cIdwcz68ReSkm6tVLtNAxuyP2sLtvAAXA4XDRwn\nn7u/oKCAW195gU0N6oHFwq/vvcF7t/63zMJvmiYTf/4CnZVMmCuAOy65keio6FLrdD67HXN7nceW\nLclER0djscjhNVE15JMmarzwyHBcOIqXnUYhkdGVnyzr1qE303tvaxquDqHNuro8fPUDJ1137h+L\n2BAbjWG1YhgG+xo34psSc+Qf9ceff3Lvy4/wY+QGtrVxsKrtER795qUy27RarcTExEjBF1VKevqi\nxht6xVBWLPybzb/vBwOS+idw8ZBLKt2uxWLhpiFjyrVuYEAAhuPYFw+mie24Yv3Mm2/z+74c8m0H\nqB/lnhraMAz2BGWSn59fxgRqQlQ9KfqixrNYLDz75lNs2rgRwzBQSUmVuvuUJ3qf153ILz9mb3gI\nOF00yHfwr+fHFz+/Z89uft9xCEtMA1ypFkyXiWFxZwzLtxMQEFCleYU4GSn6okaa/uMsNq7/h3oJ\ndbj2+iuxWCy0at26yvb/y4+zWLdsPXENYxl5y/V8M/NnMjq2ITQsFIDCnbvJyckhNNS9XFBQgMuw\nYQGCQtuy56vlRDayUS8witu6XFXlX1JCnIwUfVHjfPHpN/z81QZsRjgOVxq7d73Jg4/eWWX7//az\nb5n82A/YcgNx4mDXll1Y2ydAUcEHyIiOYPP2bcTFxQHQpElT2kcarCnMxxoYTL3cZjx7wZW0qcIv\nKiHKQ44giRply9ZtfDflVw4f2Y/LdGKzBLJpzd4qzbBi9l/Yct3j71ZsbFywmXaJTbCmHbsZXHxq\nJm3UsYnXDMNg/LgHGHN2XS5vGMSr/x4lBV/USNLTFzXGhg0befKlr7HFdSa0TiG7ty4nMbQt9sCq\n7ZvYgkpfOGUPsXNx3wtJyUhjwa5t2E2Dmy4aRmRkVKn1rFYr148YUZVRhagwKfqixpgxeyEOi/ue\ntVarnZD6zcjO2MKo66+p0hw33H09L2x+kYxNOQTEWxnxH/f+Rw8ZwegqTSKE90nRFzWGxTAwTbP4\noKfV4uT/3riLpKSqHSZJOrsVb898iy16M4mNGxETE1Ol+xfiTJIxfVFjXH/dMCLsu3EU5lGYm0qf\nrglVXvCPCg0NpX3HDmes4O/ds4OZM79id/I/Z6R9IU5GevqixoiLjeXt8eNYuOgP4uJi6dSxQ6Xa\n27BpJWu2zSc0IIaLLxxZY658XfnXTPKOPEmn9umsWhvO3n0P06VL+W+gLkRl1IzfAiGKhIaGctHA\n/pUu+CtWL2D+vkeJ6DGLglaf8d6kB72UsPL27fmEnt0yCA426NYlm4P7P63uSMKPSNEXNVZKSgrL\nV/5FdnZfuR0OAAAaV0lEQVR2hbddt3MWjTu4p00IDrNSWGe1R+2cCVaL47jlwmpKIvyRDO+IGmnK\nzJm8P/dPjgSGUc/xE8+OuQ7VokW5tzdd7v6Mo9CFxWrgyLdis5X/4+50Ovn4h+84nHuE3me3p1v7\nThV+DScTEDyIrf9spXmzAnbssmO1X+S1toU4HenpixrHNE2+nLeEgtiG2MOjOBzdkA9/mF6hNs7v\nOIolE/bj/HMbmb9vxdgdS1BQULm3v/ftl3m3IIXJIYU8sHgGv/6xsKIv46QuuPAm0rNfYt4ft3Hg\n8PP0H/gfr7UtxOlUqqevlIoHVgD9ACfwCeAC1mmtx1Y6nfBLhYWFpGfnQIkTZ/JdFbuH7fatf3DH\n6Dyio923Gly9ejXbt2uaNlUnrHvgwAFefnA8ackZxDSKYuzjt7O8MBNLYCQAuQlxzN64iv7de3n+\noo7TsdMA4uIu5+DBLK+1KUR5eNzTV0rZgAnA0bs2jwfGaa37ABal1BAv5BO13PI/f+D3Wc+xeMEk\nTNOkoKCAe6+/l9ylW3Hm5QJgyUqld5sTi/Wp5BfsJbrEfUuaNs1l794tZa47/qHX2PtbJnmbLez+\nNYO3n34Xu8NVah2bd++bLkS1qUxP/2XgXeAhwAA6aq2P/g08E+gP/FC5eKI2m/vbe3Rq/iFNO7vY\nlwK/TN9BakoiBxdk05RY9n81l8woB8F1LXxfsIe68dH06tytXG0nJvZmzZpvadfOfYfxJUsb0rNn\njzLXTdubWXxBmGEYZO0/wjUX9+STXRvIiwyjUUoGt11bubt0CVFTeFT0lVKjgBSt9a9KqaM38yz5\nV0MWEFnJbKIGcblcXj/P3SicT9NEd4+6fjwEW/6gIH8YBhYMw8AEAprUJ394AtvrRXLfmh+4PzON\n4RcMPm3bZ5/dneXLn2b27B9xugJo0+ZWoo67ZeFRsY2j2bn+sHufpklskxhuGXYVA3btYPe+vXS4\n4pziKZSF8HWe9vRHAy6lVH/gHOAzIK7E8+FAelkbHq+yN7eubrU9f05ODjePf451eelEG3bGDRrO\nAC+NbRtGSKllk2Bu/M9Ilvy0lKy/nRjRUeQ3Nwiq5+4/WJvXY+7Oddwad1W58g8efDVw9WlzvPjR\nEzx15wuk7EglvlkMT74xjoiIcOLi2nIubT17ceXky58fX84Ovp/fUx4V/aJxewCUUnOA24CXlFK9\ntdYLgEHAnPK05csHsuLiwmt9/mc+e5dFDYMxrGFkAI/++A3nnNXOK73+2Ib/4qfft9Ol3SFWb4wk\nMv4GnE47z3/xPFM+m8yshX+TedzYuqPAWZzZe++/jXGvPFK8lJ9fNZ9LX/78+HJ2qB35PeXNv9fv\nBZ5SSi0G7MBkL7YtqslhRz6G9dhUw3udeXz26SSvtH12m160Ovc71ux+g0ZtptCh88WA+6rcXgN7\ncd99o4ndW0j+umRchQ6C1u/huo4XeGXfQvirSl+cpbUu+VvYt7LtiZqlbWwCCzN3YESEYpom6Ax+\n/nMnAQGTufa6is8dn5aWyvadO2nZvDlhYeFER9chunPxH47s2rWLpx55l8MHbAQFF3D7vy4n8awE\ntu7ZSfchV9MwoaE3X54QfkeuyBWnNPrSK8id/AUfT5uL9YiNxJyW2OxhbNiQXOG2fl24kPEzfycj\nIIS4wu957KrhdGp3Tql1Jr4/mdz0uoQEAi6Y/PV8Pvvm/+jYrnJz8Qgh3OSKXHFKTqeTZhE5tHPu\nJTGvPoH2MEzTJCi44h+dT36bR05cQ+yRdUiPTWTijNknrFNY4Dpu2cThcJywnhDCM9LTFyflcrn4\nctIt9L9kEV37GEz9dhc/T+1Fg4R4/v3vuyvcXr5ZuqDnFBbyxrsTKSh00a9PV9qf05Zzz2vFto0r\nsBqROF2FnNWqDna7+6raXbuSmfjZVByFJud1acXoUXJrQiEqSoq+OKnt27eS1HYRQUHuC5eGX+kk\nOjyKy4c97VF7nRvUY8bhHIzgEIwjGezfuIU54dEYhoWlq6cy7g4XQ4dfTEhoMKtXaqLqRDP6xusA\nyMvL47Gn3yO7oAEAG7euoH79GDp17OydFyuEn5CiL04qICCIvDwr7umU3BOhWa3BHrf3wC230PD7\n79l16DCRUVFM39Qcw3APEznt9Zgzfxkd25/DgIEXMGBg6bN0tm7dwqGMIIKKdm/Y6rD0z3VS9IWo\nICn64qQSExuxdPlwkqOmUCfWyfxfmzOov+fz6BmGwchhwwDYs2c3M+d8UPycaboICLCfdNt69eoT\naM0vXnY6C6kTHVu8fCauGBaiNpLfEj/23ayveXXaE7w1+QUyMsq+gPqK4c9C3kS2r3+JoZd8R0xM\nXJnrVVSDBg25qEdTCo/sJe/IIeKD9jHm+itPun5sbCxXX9EFq2sXroJkkprk8u/bb2Dnjp3854o7\nubbrGG4fegebN232Sj4haivDNKt1+kDT16+K89X8P8+bzLKwHwlNCMB0mRyZHcxTo98snnisquze\nnUxaejpJKqn4gO2puFwuCgsLCQwMJC4unBsvu4PdczKLn4/rFsSrX718JiN7jS9/fnw5O9SK/B7/\nosrwjp/amrae0JYBABgWg9z4w2RmZhAZGVWlORo2TKRhw8Ryr2+xWAgMDCxezk7NLfV89qGc4zcR\nQpQgwzt+yu4Iwix5Y5IMO6GhYdUXyEMNW9XDaToBcJkuEpLqVnMiIWo26en7qRsuup3Xf3iC9Ij9\nkGljcItrK3QP2Zrinif/x9sh77L/nxRiG0bzn3FywzYhTkXG9CuhFowLsn37XoKDQ7CWmFStupim\nyXc//szW3ftomhDP1UOHnPIYQ214/301vy9nh1qRX8b0hWfCwmrOnOJvffwZP+o0jOAw5uzeyYFD\nH3PXzWOqO5YQtYoUfVFjrPxnD0ZwPACWwFD+2r63XNst/us3Fm79EtOWT4zrbG4eMa7Kz0ISwlfI\ngVxRYSkp+/hp1kcsXDwdbw4PBlhKF+rAclxslZWVyfzkCTS8MJPEPvlYOi1n2uxPvJZJiNpGir6o\nkB07NzN9yWga93wHV4PH+GLyw15re/RlAwjLTCY/4zAhmcmMurTfabfZvSeZ4IZHipeDIqyk5+/z\nWiYhahsZ3hEVsmzNZ3TudwgwqBMH++N/4+DBu4mLq/yVuud17sRnSS3ZvmMHTRo3Jjw84rTbNG7U\nlJyfoqGJe4qGzP0OzopSlc4iRG0lRd/POZ1O/u/r99mcd5gI08Y9g6+lcYNGp9ii9HBOgTOPeyff\nT8Og+tw3/D6iIqMrlScsLJy2bcp/M/KQkBCGd3iA2XM/wbQV0CCoPQMHXV6pDELUZlL0a4lFi5ey\ndMkagkPs3DjmGoKCgsq13RuTP+XH2CwsQe6zeB6Z8j6f//eZk65/bpuRLJi7nA59D5OR6mLhtmjC\n+8eSbhbw9k/v8PBI7w33lFerlufQquWrVb5fIXyRFP1aYO68hbw7YQ6GJRaXK4/NW17gtfGPl+sM\nlh25aVjqHpvWYI+1gIKCAgICAspcv2nTVgQFf8SyBdP59u+fiboyCXDPoJltO1LmNkKImkMO5NYC\nfyxeg2FxTzNssVjZtcvBwYMHy7VtvC0Ys8TtCOMctpMW/KPq10tkyODbaJTQrvgxZ76DOGI8SC+E\nqEpS9GuIP//+i3e+/oJ5SxZXeFub3Sh16qTN6iQkJKRc294zYgw9tzupuyGFlhvSeXjgyHLv976h\n91J/WRRhK+y0XJfIHZf/t8LZhRBVS4Z3aoDvZk3n7fVryY+JwfrnH/xrz25uHnFVube/6car2PrQ\neA6kBGK35zN0SEfCwso3eVpQUBAv3Xy/R7nrRNfhkese8WhbIUT1kKJfA8zYsI78GPfQiDMqil+2\naG6uwPYxMTG8/dbjbNu2jdjY2DJPn3zli4nM37MLmwHXtO/CFQMGeyl91XE63bNp1oR5goTwVR4V\nfaWUDZgINAECgGeBDcAnuG+ouk5rLdMdlpPluKtaPRlzCwgIoFWrVmU+98Nvv/BtThpm4wQA3ly7\ngk4tWxEX167M9U9nx/aNbFr3Pjabg9i6l9K+4wCP2qmIl75+l3mZm7GYBhfFt2Ps8FFnfJ9C1Eae\njumPBA5prXsDg4C3gPHAOK11H8CilBripYy13hXndiVk/35Ml4uAgwcZfk5Hr7a/LWU/ZsSxidVy\nYqNZu2WTR21lpKexbd2dXDloDsP7LyDQ+Th645/eilqmGQt+5ZfoA+S3r09uh3pMs2xhycplZ3Sf\nQtRWnhb9b4FHS7ThADpqrRcWPTYTOP019AKAi3qfz7tXXsvY6BjevmQoVw++xKvtd26uCDiUWrwc\nfeAQXdt18KitdWsXcmHPY9McdO2Qw64d8yud8VSSD++D6GMHps26EWzbu+OM7lOI2sqjoq+1ztFa\nH1FKhQPfAQ8DJU8KzwIivZDPb7RsdhbXDxtBm1atvd527y7n8d/mbWh3MIMOBzN4ot8l1Ktbz6O2\n6ic0Z+v2Y/eyTU03sQee2btV9Wp7LoFbjn1phW1MpU/H7md0n0LUVh7fREUplQhMBd7SWn+qlNql\ntW5U9NxlQD+t9enO4avWO7gIz8z46XVy0j4lMMBBek5fRo5+/YxPZfz7koV8s3IuFmBUj8Gc177z\nGd2fEDWcx79wHhV9pVRdYC4wVms9t+ixH4BXtNYLlFLvAnO01t+dpim5c1Y1qkx+p9OJw+EodZPy\nqubP73918+XsUCvyV/mdsx4CooBHlVKP4e6x3wm8qZSyAxuByZ6GEpXz89QZzP1pCYYBg6++kAsG\n9PX6PqxWq5w6KYQP8qjoa63vAu4q46m+lUojKu2v5X8xafyvWHLdF2d9/M80GjdL5KzmZ1VzMiFE\nTSDTMNQyq1esLS74AGSG8+cfcnqjEMJNin4tc1bLZjhtOcXLrqBs2pzTxuP2tm3fwPcz3mHhH9O9\nEU8IUc2k6NcyHbu0p8GAdJyt/oKEFIb8uyftzin/TUlK+nvNQpbuuIMGvT8np+6zTPrh+Qptv3fv\nHlatWkV+fr5H+xdCeJ/MvVOL5Obm8vbUsXS76xAWi8Hm3w7Tp7fn57Ov3z4ZdaH7r4Y6dU12rP8N\nh+M+bLbTf2w+/HwSU+dvwmEEkxAyhZcevZO42FiPswghvEN6+rXIomWzaXphChaL+2yuFhceYd6K\naR63d8LJvCblOh8/KyuTH+avwxregMCwOhwyEvnwi9OdvSuEqApS9GuR0KBw8o64ipedhSZrVm3k\n28nfU1hYWOH2OrS4lvWLwzFNk/07DeJsl5TrNM2cnFwc5rGrdg3DoNBxig2EEFVGin4t0q3L+WQt\n78Ch5HwyDxXw4yt5rNvaki9n7ubuB5/D4ahY5W3Tuiv9237A/kW3kFDwAsMvvrtc28XHx9Mi3orL\n6d6fJS+Fvt3bV/j1CCG8z+NpGLxErsj1MtM0WbV2BT/+NJ2//4nFYnX3uB2Fedw0ogWXXXxsHv0z\nmT8/P5+JX35Ddk4Bvbq257xzvT9tQk18/yvCl/P7cnaoFfmr/IpcUUMZhkGHdueycuUWVu3IKP1k\nFX6/BwYGcvuYG6puh0KIcpHhnVrqyssvIcq+G9N04XI6SIxOY9DA/tUdSwhRzaSnX0uFh0fw5osP\nMe2H6djtNoYPvRm73X76DYUQtZoUfR+yePESdmzfRY+e59GkSePTrh8WFsb115X/ButCiNpPir6P\neG/CZ/w2exdWI4IZP3/IXf8byrldOp2wnsvl4rOfPySlcB+RRh3GXHqb9PCFEMVkTN8HuFwu5s1Z\nj80SiWEYuBzxfP/93DLXnTD1DTY1X07WuQfY0XYdr337QhWnFULUZFL0fZRpln3G1l7XTuxh7p69\nLdBKim1vVcYSQtRwUvR9gMVioVuP5jid2QAY1kMMvrjsOXUCHEGllgMLg8pcTwjhn2RM30f8986b\nad16DsnJ++nWrT9JrZLKXG9kj5t45/dXyIvOIiA9hNHdRlVtUCFEjSZF34f063/Baddp1vgsXhr1\nNllZmYSHR5zxG5YLIXyLFP1ayDAMIiIiqzuGEKIGkqJfi2i9mW8/+wmXEy4c3I2elZhLXwhRO0nR\n9wHL16xg7T+bOKf52XRq06HMdVJTU3n+ofdxprtvVLJ59Q9ERITTrr1nd80SQtROcvZODTfp12k8\ntW0KUxKTeXLzN0yZ81OZ6y1ZvJT81IjiZUt+NH8sXF5VMYUQPkKKfg33S/IKXI3cxdzZOJKZO5eV\nuV7jJo1wWY8ULztcecTGRVdJRiGE75CiX8Mdf+7Nyc7FadO2Df1HnI0jcD+Ftv2c3SuUy68aeqbj\nCSF8jFfH9JVSBvAOcA6QB9yktf7Hm/vwN4MadeHTXUswG0Vi3ZHJxU17nXTdm26/gX/ddA0Oh4Pg\n4OAqTCmE8BXePpA7FAjUWndXSnUFxhc9Jjx0Rb8hNFvXiDXbNtCxVVvOad3ulOvb7XaZYE0IcVLe\nLvo9gVkAWus/lVLev0deDbJ561ZWrdGc26EDYWHhZ2w/ndp0OOlZO0IIURHeLvoRQMl79DmUUhat\ntcvL+6l2H036lkkrNlMYEEb8D7/x4h1jaNK4SXXHEkKIU/J20c8ESnZ5T1vw4+LOXA/5TCkoKOD7\nZWsxoxtjA1KDw/hq5i+8+ui91R2twnzx/S9J8lcfX84Ovp/fU94u+ouBS4DJSqnzgLWn28AX70h/\n5MgRClylz6PJznX43GuJiwv3ucwlSf7q48vZoXbk95S3T9mcBuQrpRYDrwB3e7n9GiE0NJT29SJw\nFuQBYM86QL9z21dzKiGEOD2v9vS11iZwuzfbrKmeu+9uvp81nV37U+ne/iK6djrx1oVCCFHTyNw7\nHrJardw26jqf/hNRCOF/5IpcIYTwI1L0hRDCj0jRF0IIPyJFvxqYpsnGTZtYuepvnE5ndccRQvgR\nOZBbxUzT5NG33uT39GxMq412P//EW/c/SEBAQHVHE0L4AenpV7Ely5fx+5F8rLHx2KLrsD4yns+m\nTanuWEIIPyFFv4qlZWZA4LFpjw2bjdzCwmpMJITwJ1L0q9gFPXrRJO0Apss9JVH0/mQu7d2nmlMJ\nIfyFFP0qFhwczIR77uHyUDtDgiy8NnoUTRo1ru5YQgg/IQdyq0FkZBT/GzW6umMIIfyQ9PSFEMKP\nSNEXQgg/IkVfCCH8iBR9IYTwI1L0hRDCj0jRF0IIPyJFXwgh/IgUfSGE8CNS9IUQwo9I0RdCCD8i\nRV8IIfyIFH0hhPAjHk24ppSKAL4AIgA78D+t9Z9KqfOA14BC4Fet9VNeSyqEEKLSPO3p/w/4TWvd\nFxgNvFP0+LvA1VrrXkBXpVT7ykcUQgjhLZ5OrTweyC/62Q7kKqXCgQCt9Y6ix38BLgRWVSqhEEII\nrzlt0VdKjQHuBkzAKPp3tNb6L6VUPeBz4L+4h3oyS2yaBTT1emIhhBAeO23R11pPBCYe/7hSqi3w\nFXCP1npRUU8/osQq4UC6t4IKIYSoPMM0zQpvpJRqDUwBrtRary3x+ErgcmAH8DPwhNZ6uXeiCiGE\nqCxPx/SfAwKB15VSBpCutR4G3I67928BZkvBF0KImsWjnr4QQgjfJBdnCSGEH5GiL4QQfkSKvhBC\n+BEp+kII4Uc8PXvHY7Vh3p6iM5beAc4B8oCbtNb/VG+qU1NK2XBfb9EECACeBTYAnwAuYJ3Wemx1\n5SsPpVQ8sALoBzjxrewPApfh/sy/AyzAR/IXfXY+xf3ZcQA34yPvv1KqK/CC1vp8pdRZlJFZKfUY\ncDHu2nN3TTrr8Lj87YE3cP8f5AM3aK0PKqVuBm7Bnf9ZrfX0U7VZHT392jBvz1AgUGvdHXgI97QU\nNd1I4JDWujcwCHgLd+5xWus+gEUpNaQ6A55KUeGZAOQUPeRL2fsA3Yo+L32BRvhQfmAwYNVa9wCe\nxn3Kdo3Pr5S6D/gA9+nlUEZmpVQHoLfWuitwDfB29aQ9URn5XwPGaq0vAKYBDyil6gJ3AN2Ai4Dn\nlVL2U7VbHUV/PPBe0c+nm7enpuoJzALQWv8JdK7eOOXyLfBo0c8W3L2FjlrrhUWPzcTdg66pXsbd\nMdiLezoQX8o+EFinlPoe+BH3hYu+lH8zYCv6CzcSd4/SF/JvBYaVWO50XOb+uH+XZwNorZMBq1Iq\npkpTntzx+a8qcTGsDfcoQxdgkdbaobXOBLYA7U7V6Bkt+kqpMUqptUqpNUf/BVporfNLzNvzIGXP\n2xN5JrNVUgSQUWLZoZSq0cdHtNY5WusjRV+w3wEP4y6eR9XY91wpNQpI0Vr/yrHMJd/vGpu9SCzQ\nCRiB+wLGL/Gt/Nm459HahLvD9gY+8NnRWk/D3bk5qqzM4ZT+Xc6mhryW4/NrrQ8AKKW6A2OBVzmx\nFp02/xkd06/F8/Zk4s54lEVr7aquMOWllEoEpgJvaa0nKaVeLPF0TX7PRwMupVR/3MdRPgPiSjxf\nk7MDHAY2aq0dwGalVB7QsMTzNT3/3cAsrfXDSqkGwDzcx4WOqun5jyr5OxoOpOH+XfaZ2qOUugr3\nkPJgrfVhpVSF81d577Ro3p5vgWu11kf/rMoC8pVSTYv+hBwILDxFM9VtMe5xTooOQK899erVr2js\n7xfgfq31p0UP/62U6l308yBq6Huute6jtT5fa30+7qm6rwdm+kL2Iotwj7eilEoAQoHfi8b6oebn\nT+VYbzIdd2fxbx/Kf9TKMj4zfwADlFKGUqoRYGitU6st4SkopUbi7uH31VrvLHp4GdBTKRWglIoE\nkoB1p2qnys/eoXbM2zMN6K+UWly0PLo6w5TTQ0AU8GjR2QomcCfwZtGBn43A5GrMV1H3Ah/4Qnat\n9XSlVC+l1DLcQwy3456U8ENfyI/7AOJEpdQC3MfhHgT+wnfyH3XCZ0ZrbSqlFgJLcP/f1NSzkCzA\n68BOYJpSygTma62fVEq9gbtjYeA+UF1wqrZk7h0hhPAjNfrgoxBCCO+Soi+EEH5Eir4QQvgRKfpC\nCOFHpOgLIYQfkaIvhBB+RIq+EEL4ESn6QgjhR/4fRtN7tNTbaIIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x114ef5160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(c[:,0],c[:,1],c=c_space,cmap=cmap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        ,  0.94772484],\n",
       "       [ 0.94772484,  1.        ]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# the correlation coefficient is still relatively large\n",
    "np.corrcoef(c,rowvar=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try relating with the nosie scale and sparsity of the sequences with the correlation coefficient.\n",
    "\n",
    "\n",
    "# Ordinal comparison\n",
    "\n",
    "e.g. comparing rank"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Get Tweets from the Twitter public API\n",
    "\"\"\"\n",
    "import yaml\n",
    "\n",
    "# Get your Twitter API tokens\n",
    "# this is specific to my computer; modify for yours\n",
    "my_creds_file = '/Users/jkolb/.twitter_api_creds'\n",
    "creds = yaml.load(open(my_creds_file))\n",
    "consumer_key = creds['audience']['consumer_key']\n",
    "consumer_secret = creds['audience']['consumer_secret']\n",
    "access_token_key = creds['audience']['token']\n",
    "access_token_secret = creds['audience']['token_secret']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "api = twitter.Api(consumer_key=consumer_key,\n",
    "                 consumer_secret=consumer_secret,\n",
    "                 access_token_key=access_token_key,\n",
    "                 access_token_secret=access_token_secret\n",
    "                 )\n",
    "\n",
    "mom_tweets = []\n",
    "for _ in range(20):\n",
    "    mom_tweets.extend( api.GetSearch(\"mom\",count=100) )\n",
    "    time.sleep(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Get Tweets from the Gnip Search API\n",
    "\"\"\"\n",
    "\n",
    "\"\"\"\n",
    "from search.api import Query\n",
    "import json\n",
    "import yaml\n",
    "creds = yaml.load(open('/Users/jkolb/.creds.yaml'))\n",
    "\n",
    "# set up a query to the Gnip Search API\n",
    "q = Query(creds['username'],\n",
    "          creds['password'],\n",
    "          creds['search_endpoint'],\n",
    "          paged=True,\n",
    "          hard_max = 1000,\n",
    "          search_v2 = True\n",
    "          )\n",
    "\n",
    "# query parameters\n",
    "start_date = '2016-08-01T00:00'\n",
    "end_date = '2016-09-01T00:00'\n",
    "rule = 'mom'\n",
    "\n",
    "# get the tweet data\n",
    "q.execute(rule,start=start_date,end=end_date)\n",
    "mom_tweets = list(q.get_activity_set())\n",
    "\"\"\"\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def token_frequency(tweets,num_top_terms = 30,delta=10**-5,epsilon=0.001,):\n",
    "    \"\"\"\n",
    "    Space-tokenize tweet bodies and \n",
    "    return exact and approximate 1-gram counts\n",
    "    \n",
    "    Approximate the counts with a count-min sketch\n",
    "    \"\"\"\n",
    "    \n",
    "    terms = collections.defaultdict(int)\n",
    "    sketch = count_min.Sketch(dict(delta=delta,epsilon=epsilon,k=num_top_terms))\n",
    "    for tweet in tweets:\n",
    "        try:\n",
    "            tokens = tweet['body'].split()\n",
    "        except TypeError:\n",
    "            tokens = tweet.AsDict()['text'].split()\n",
    "        for token in tokens:\n",
    "            terms[token.lower()] += 1\n",
    "            sketch.update(token.lower(),1)\n",
    "    return (terms,sketch)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def display_two_lists(list_1,list_2):\n",
    "    \"\"\"helper function\"\"\"\n",
    "    for x,y in zip(list_1,list_2):\n",
    "        print(x,y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "num_top_terms = 30\n",
    "# accuracy parameters for CM sketch\n",
    "delta = 10**-4\n",
    "epsilon = 0.01"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# get exact and approximate top terms and counts\n",
    "mom_terms,mom_sketch = token_frequency(mom_tweets,num_top_terms,delta,epsilon)\n",
    "exact_top_mom_terms = list(reversed(sorted(mom_terms.items(),key = operator.itemgetter(1))))[:num_top_terms]\n",
    "approx_top_mom_terms = [(term,count) for count,term in reversed(sorted(mom_sketch.top_k.values(),key = operator.itemgetter(0)))]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Kendall's tau coefficient** is a sort of correlation coefficient that is proportional to the difference between the number of _concordant_ pairs and the number of _discordant_ pairs.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/jkolb/pyenvs/3.5.1/lib/python3.5/site-packages/scipy/stats/stats.py:250: RuntimeWarning: The input array could not be properly checked for nan values. nan values will be ignored.\n",
      "  \"values. nan values will be ignored.\", RuntimeWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.78551229330666639"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "kt_result = scipy.stats.kendalltau(exact_top_mom_terms,approx_top_mom_terms)\n",
    "kt_result.correlation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('mom', 1507) ('mom', 1512)\n",
      "('rt', 1236) ('rt', 1248)\n",
      "('my', 941) ('my', 946)\n",
      "('to', 901) ('to', 911)\n",
      "('and', 745) ('and', 783)\n",
      "('your', 710) ('your', 720)\n",
      "('i', 648) ('i', 648)\n",
      "('you', 551) ('you', 551)\n",
      "('a', 431) ('mom:', 453)\n",
      "('mom:', 424) ('a', 431)\n",
      "('me', 415) ('the', 425)\n",
      "('the', 413) ('me', 419)\n",
      "('this', 283) ('is', 299)\n",
      "('is', 279) ('this', 286)\n",
      "('me:', 231) ('me:', 247)\n",
      "('when', 229) ('when', 236)\n",
      "('how', 219) ('her', 219)\n",
      "('her', 205) ('how', 219)\n",
      "('for', 201) ('for', 201)\n",
      "('of', 185) ('of', 185)\n",
      "('she', 175) ('she', 183)\n",
      "('it', 160) ('been', 178)\n",
      "('been', 157) ('through', 177)\n",
      "(':', 156) ('@pearcegodeep:', 172)\n",
      "('many', 156) ('\"', 167)\n",
      "('\"', 156) ('it', 164)\n",
      "('through', 150) ('https://t.co/esugkx28mv', 164)\n",
      "('@pearcegodeep:', 142) (':', 160)\n",
      "('https://t.co/esugkx28mv', 142) ('many', 156)\n",
      "('do', 132) ('do', 154)\n"
     ]
    }
   ],
   "source": [
    "display_two_lists(exact_top_mom_terms,approx_top_mom_terms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A harder problem: how to account for the varying importance of rank?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Over/under-indexing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# get some data around the term 'dad'\n",
    "\n",
    "# public API version\n",
    "dad_tweets = []\n",
    "for _ in range(20):\n",
    "    dad_tweets.extend( api.GetSearch(\"dad\",count=100) )\n",
    "    time.sleep(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"\\n# Gnip API verions\\nrule = 'dad'\\nq.execute(rule,start=start_date,end=end_date)\\ndad_tweets = list(q.get_activity_set())\\n\""
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "# Gnip API verions\n",
    "rule = 'dad'\n",
    "q.execute(rule,start=start_date,end=end_date)\n",
    "dad_tweets = list(q.get_activity_set())\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dad_terms,dad_sketch = token_frequency(dad_tweets,num_top_terms,delta,epsilon)\n",
    "exact_top_dad_terms = list(reversed(sorted(dad_terms.items(),key = operator.itemgetter(1))))[:num_top_terms]\n",
    "approx_top_dad_terms = [(term,count) for count,term in reversed(sorted(dad_sketch.top_k.values(),key = operator.itemgetter(0)))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('dad', 1549) ('mom', 1507)\n",
      "('rt', 1127) ('rt', 1236)\n",
      "('to', 980) ('my', 941)\n",
      "('i', 860) ('to', 901)\n",
      "('the', 844) ('and', 745)\n",
      "('me', 714) ('your', 710)\n",
      "('you', 690) ('i', 648)\n",
      "('your', 675) ('you', 551)\n",
      "('on', 675) ('a', 431)\n",
      "('he', 655) ('mom:', 424)\n",
      "('my', 655) ('me', 415)\n",
      "('before', 538) ('the', 413)\n",
      "('please', 526) ('this', 283)\n",
      "('yo!', 510) ('is', 279)\n",
      "('planet.', 510) ('me:', 231)\n",
      "('@donaldjtrumpjr!', 510) ('when', 229)\n",
      "('destroys', 510) ('how', 219)\n",
      "('ask', 510) ('her', 205)\n",
      "('twitter.', 510) ('for', 201)\n",
      "('follow', 510) ('of', 185)\n",
      "('resign', 510) ('she', 175)\n",
      "('noticed', 510) ('it', 160)\n",
      "('thanks…', 490) ('been', 157)\n",
      "('@sethrogen:', 490) (':', 156)\n",
      "('and', 376) ('many', 156)\n",
      "('is', 346) ('\"', 156)\n",
      "('a', 339) ('through', 150)\n",
      "('for', 231) ('@pearcegodeep:', 142)\n",
      "('are', 193) ('https://t.co/esugkx28mv', 142)\n",
      "('this', 182) ('do', 132)\n"
     ]
    }
   ],
   "source": [
    "display_two_lists(exact_top_dad_terms,exact_top_mom_terms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But we don't really care about rank here. We care about removing the effect of a baseline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def normalize_by_difference(term_counts,baseline_counts):\n",
    "    \"\"\"define a normalized term frequency that subtracts off a baseline count\"\"\"\n",
    "    normed_term_counts = {}\n",
    "    for term,count in term_counts.items():\n",
    "        try:\n",
    "            normed_term_counts[term] = count - baseline_counts[term]\n",
    "        except KeyError:\n",
    "            normed_term_counts[term] = count\n",
    "    return normed_term_counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('dad', 1462),\n",
       " ('on', 590),\n",
       " ('he', 563),\n",
       " ('before', 538),\n",
       " ('yo!', 510),\n",
       " ('planet.', 510),\n",
       " ('resign', 510),\n",
       " ('@donaldjtrumpjr!', 510),\n",
       " ('destroys', 510),\n",
       " ('twitter.', 510),\n",
       " ('noticed', 510),\n",
       " ('ask', 494),\n",
       " ('please', 492),\n",
       " ('thanks…', 490),\n",
       " ('@sethrogen:', 490),\n",
       " ('follow', 484),\n",
       " ('the', 431),\n",
       " ('me', 299),\n",
       " ('i', 212),\n",
       " ('you', 139)]"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "normalized_results = normalize_by_difference(dad_terms,mom_terms)\n",
    "\n",
    "# look at top of list to see most \"dad\"-like terms\n",
    "list(reversed(sorted(normalized_results.items(), key = operator.itemgetter(1))))[:20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('go', -77),\n",
       " ('at', -77),\n",
       " ('did', -86),\n",
       " ('a', -92),\n",
       " ('been', -95),\n",
       " ('times', -97),\n",
       " ('do', -99),\n",
       " ('this', -101),\n",
       " ('an', -102),\n",
       " ('someone', -105),\n",
       " ('rt', -109),\n",
       " ('her', -114),\n",
       " ('many', -153),\n",
       " ('\"', -155),\n",
       " ('when', -160),\n",
       " ('how', -191),\n",
       " ('me:', -227),\n",
       " ('my', -286),\n",
       " ('and', -369),\n",
       " ('mom', -1412)]"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# and look at the bottom of the list to see the mom-like terms\n",
    "list(reversed(sorted(normalized_results.items(), key = operator.itemgetter(1))))[-20:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But this normalization-by-difference only works if the two term frequency distributions have the same scale of counts. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def normalize_by_fraction(term_counts,baseline_counts):\n",
    "    \"\"\"normalize counts by the difference of the term-fractions for each distribution\"\"\"\n",
    "    normed_term_counts = {}\n",
    "    B = sum(baseline_counts.values())\n",
    "    A = sum(term_counts.values())\n",
    "    for term,count in term_counts.items():\n",
    "        try:\n",
    "            # fraction of baseline distribution for this term\n",
    "            b_frac = baseline_counts[term]/B\n",
    "        except KeyError:\n",
    "            b_frac = 0\n",
    "        \n",
    "        # fraction of primary term frequency distribution for this term\n",
    "        a_frac = count/A\n",
    "        \n",
    "        factor = (a_frac-b_frac)\n",
    "        normed_term_counts[term] = count * factor \n",
    "    return normed_term_counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('dad', 65.04625878816911),\n",
       " ('on', 11.35391535337112),\n",
       " ('he', 10.49533948411985),\n",
       " ('the', 9.62448294280985),\n",
       " ('before', 8.356979933593186),\n",
       " ('yo!', 7.50974447812906),\n",
       " ('planet.', 7.50974447812906),\n",
       " ('resign', 7.50974447812906),\n",
       " ('@donaldjtrumpjr!', 7.50974447812906),\n",
       " ('destroys', 7.50974447812906),\n",
       " ('twitter.', 7.50974447812906),\n",
       " ('noticed', 7.50974447812906),\n",
       " ('please', 7.42691656924851),\n",
       " ('ask', 7.25358374982895),\n",
       " ('follow', 7.0934832946413815),\n",
       " ('thanks…', 6.932293922332899),\n",
       " ('@sethrogen:', 6.932293922332899),\n",
       " ('me', 5.417262730735218),\n",
       " ('i', 3.85985063825802),\n",
       " ('you', 1.8111924301860725)]"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "normalized_results = normalize_by_fraction(dad_terms,mom_terms)\n",
    "sorted_list = list(reversed(sorted(normalized_results.items(), key = operator.itemgetter(1))))\n",
    "sorted_list[:20]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Well, this sort of works. We could also comine the term-fractions in other ways, e.g. a ratio (a_frac/b_frac), or the relative difference ( [a_frac-b_frac]/a_frac ). \n",
    "\n",
    "In the end, we need to think harder about what differences and similarities we want to be able to highlight and ignore. See Part Two..."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
